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ABSTRACT 


For  the  spacecraft  and  aircraft  designers,  the  ability  to  change  and  control  the  shape  of 
the  structures  has  been  a  challenging  problem.  For  a  spacecraft,  it  is  highly  desirable  to 
change  the  reflector  shape  in  orbit  to  conqwnsate  for  the  sur&ce  errors  and  also  to 
perform  antenna  beam  shaping.  In  the  case  of  aircraft,  the  sh£^  control  of  propellers, 
helicopter  rotor  blades,  and  aircraft  wings  can  increase  efficiency  and  maneuverability 
and  reduce  vibration  and  noise  of  the  vehicles.  In  the  present  work,  the  sht^  control  of 
fiber  reinforced  composite  plates  with  embedded  piezoelectric  actuators  is  investigated. 

In  the  present  study,  a  finite  element  formulation  is  developed  for  modeling  a 
laminated  composite  plate  that  has  distributed  piezoelectric  actuators  and  sensors 
subjected  to  both  mechanical  and  electrical  loads.  A  sinq>le,  h^er  order,  shear 
deformation  theory  with  Hamilton's  principle  is  used  to  formulate  the  equations  of 
motion.  The  model  represents  the  parabolic  distribution  of  transverse  shear  stresses  and 
the  non-linearity  of  in-plane  displacements  across  the  thickness.  The  model  is  valid  for 
both  segmented  and  continuous  piezoelectric  elements,  which  can  be  either  surface 
bonded  or  embedded  in  the  laminated  plate.  A  four-node,  bilinear,  isoparametric, 
rectangular  element  with  seven  degrees  of  fi'eedom  at  each  node  is  developed.  The 
electric  potential  is  treated  as  a  generalized  electric  coordinate  like  the  generalized 
di^lacement  coordinates  at  the  mid-plane  of  the  actuator  layers. 

For  shape  control,  an  optimization  algorithm,  based  on  a  finite  element  techniques,  is 
presented  for  an  optimal  jq)plied  voltage  to  each  actuator  to  minimize  the  error  between 
the  desired  shape  and  the  actual  shape.  The  error  (objective)  function  is  the  mean  square 
of  the  error  between  the  point  in  the  actual  surftice  and  the  corresponding  point  in  the 
desired  surfece.  Based  on  these  techniques,  two  computer  programs  have  been  developed, 
a  finite  element  modeling  of  a  composite  plate  with  piezoelectric  actuators  and  an 
optimization  model  of  the  actuator  voltages  for  shape  control.  The  present  work 
demonstrates  the  feasibility  of  the  application  of  the  piezoelectric  actuators  for  the  shape 
control  of  con^site  plates  used  in  aero^ace  structures. 
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I.  INTRODUCTION 


Advanced  structures  with  integrated  self-monitoring  and  control  capabilities  are 
becoming  increasingly  iitqx)rtant  due  to  the  rapid  development  of  ‘intelligent’  space 
structures.  Those  smart  structures  that  have  actuators  distributed  throughout  are  defined 
as  adaptive,  or  actuated,  structures.  Exan5)les  of  such  adaptive  structures  are  conventional 
aircraft  wings  with  articulated  leading-and  trailing-edge  control  sur&ces,  robotic  systems 
with  articulated  manipulators  and  end  effectors  and  spacecraft  antenna  reflectors. 
Structures  that  have  sensors  distributed  throughout  are  a  subset  referred  to  as  sensory. 
These  structures  have  sensors  which  can  detect  displacements,  strains  or  other  mechanical 
states  or  properties  like  electromagnetic  states  or  properties,  ten^)erature,  heat  flow,  or 
damage.  Applications  of  this  technology  might  include  damage  detection  in  long  life 
structures,  or  embedded  or  conformal  RF  antennas  within  a  structure.  The  overlap 
structures,  winch  contain  both  actuators  and  sensors  implicitly  linked  closed-loop 
control,  are  referred  to  as  controlled  structures.  A  subset  of  controlled  structures  are 
active  structures,  distinguished  fi“om  controlled  structures  by  extensively  distributed 
actuators,  which  have  structural  functionality  and  are  part  of  the  load  bearing  system. 

Intelligent  structures  are  a  subset  of  active  structures  that  have  highly  distributed 
actuator  and  sensor  systems  with  structural  functionality  and,  in  addition,  distributed 
control  functions  and  computing  architecture  [Ref  1]. 

Smart  structures  have  applications  due  to  their  ability  to  change  shape.  Some  of  these 
are  for  spacecraft  antennas,  to  compensate  for  surfece  error  and  thermal  distortion  to 
inprove  antenna  performance,  and  to  change  the  antenna  beam  shape  in  orbit.  They  may 
be  used  for  and  submarine  and  helicopter  shape  control  as  weU.  An  aeroelastic  application 
to  aircraft  structures  is  quasi-static  control  of  camber,  dynamic  control,  and  flutter 
suppression.  Smart  structures  can  be  used  in  acoustic  control  by  developing  an  adaptive 
structure  in  which  the  structural  response  can  be  modified  with  var5Tng  input  disturbances. 
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The  piezoelectric  material  is  one  of  the  smart  materials  that  can  be  used  as  an 
actuator  and  sensor.  Piezoelectricity  refers  to  the  phenomenon  of  generating  an  electric 
charge  in  a  material  when  subjecting  it  to  mechanical  stress,  which  is  known  as  a  direct 
piezoelectric  effect.  The  converse  piezoelectric  effect  is  described  as  an  induced  strain  in 
response  to  an  applied  electric  field. 

Piezoelectric  properties  occur  naturally  in  some  crystalline  materials  and  can  be 
induced  in  other  polycrystalline  materials.  The  distortion  of  the  crystal  domain  produces 
the  piezoelectric  effect.  The  domain  may  be  aligned/poled  by  the  q>plication  of  a  large 
electric  field,  usualfy  at  high  tenqjerature.  Subsequent  ^plication  of  the  electric  field  will 
produce  additive  strains  in  the  local  domain,  >\iiich  translate  into  a  global  strain  in  the 
material. 

The  piezoelectric  effect  was  discovered  in  1880  by  Pierre  and  Jacques  Curie,  The 
direct  piezoelectric  effect  has  been  used  for  a  long  time  in  sensors  such  as  accelerometers. 
Use  of  the  converse  effect,  however,  has  imtfl  recently  been  restricted  to  ultrasonic 
transducers.  Barium  titanate,  discovered  in  1940,  was  the  first  widely  used  piezoceramic. 
Lead  zirconate  titanate  (PZT)  [Ref  2],  discovered  ini 954,  has  now  largely  superseded 
barium  titanate  because  of  its  stronger  piezoelectric  effect. 

A  P-phase  polymeric  piezoelectric,  polyvinylidene  fluoride  (PVDF),  was  initially 
discovered  by  Kawai  in  1969  [Ref  3].  Raw  polymeric  PVDF  (a-phase)  is  an  electrical 
insulator,  and  it  does  not  have  any  intrinsic  piezoelectric  properties.  If  the  raw  material 
is  polarized  during  the  manufecturing  process,  PVDF  transforms  to  a  P-phase  ,  a  tough 
and  flexible  semi-crystalline  material,  and  can  be  made  to  strain  either  in  one  or  two 
directions  in  the  film  plane.  Since  P-phase  PVDF  possesses  a  strong  direct  piezoelectric 
effect,  it  has  been  used  in  many  transducer  applications;  e.g.,  sonar,  medical  ultrasonic 
equipment,  robot  tractile  sensors,  acoustic  pick-ups,  force  and  strains  gages,  etc.  Due  to 
its  distinct  characteristics,  such  as  flexibility,  durability,  manufecturability,  etc,,  PVDF  is 
commonly  used  in  such  structural  systems. 
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As  a  first  step  in  the  following  work,  a  finite  element  analysis  of  a  graphite/epoxy, 
fiber  reinforced,  con^site,  laminated  plate  is  developed,  using  a  single  higher-order 
shear  deformation  theory  [Ref  4].  Multi-layered  con^wshes  have  found  wide  use  in  many 
weight-sensitive  structures,  such  as  aircraft,  spacecraft  antennas,  and  missile  structural 
conqwnents,  where  high  strength-to-weight  and  stifl5iess-to-weight  ratios  are  required,  A 
laminate  is  a  multi-layered  conyx>site  made  up  of  several  individual  layers  (laminae)  in 
vvluch  the  fibers  are  oriented  in  a  predetermined  direction  to  provide  eflBciently  the 
required  strength  and  stifibess  parameters  in  each  laminae. 

For  analysis  and  design  of  structural  laminates,  a  classical  plate  theory  (CPT)  [Refe. 
5,6,7,  &  8]  has  been  used  in  t^ch  it  is  assumed  that  normals  to  the  mid-plane  before 
deformation  remain  straight  and  normal  to  the  mid-plane  after  deformation  (classical 
Kirchhoff  hypotheses).  This  theory  under-predicts  deflections  and  over-predicts  natural 
fi'equencies  and  bucklii^  loads.  These  results  are  due  to  the  neglect  of  transverse  shear 
strains  in  the  classical  theory.  These  errors  are  even  higher  for  a  plate  made  of  advanced 
conqKtshes  like  graphite  epoxy  and  boron-epoxy,  whose  elastic  modulus  to  shear 
modulus  ratio  are  very  large  (Le.,  of  the  order  of  25  to  40  instead  of  2.6  for  a  typical 
isotropic  materials).  These  High  ratios  render  classical  theories  inadequate  for  the  analysis 
of  conqwshe  plates.  An  adequate  theory  must  account  for  transverse  shear  strains. 

Many  plate  theories  have  been  proposed  to  incorporate  the  influence  of  the  transverse 
shear  strains.  In  one  of  these,  the  Reissner-Mindlin  plate  theory  [Ref  9],  transverse  shear 
and  rotary  inertia  effects  are  included,  and  it  contains  a  displacement  field  which  accounts 
for  linear,  or  higher-order,  variations  of  mid-plane  displacements  through  the  thickness, 
but  on  the  other  side  the  deviation  increases  with  increasing  mode  numbers.  A  theory  for 
laminated  isotropic  plates  by  Stavsky  [Ref  10].  has  been  generalized  to  laminated 
anisotropic  plates  by  Yang,  Norris  et.  al.  [Ref  11].  Whitney  [Ref  12]  has  presented  an 
approximate  method  to  incorporate  the  influence  of  shear  deformation  on  plate  deflection, 
in  a  flexural  vibration  and  buckling  analysis.  Elasticity  solutions  by  Pagano  and  his 
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associates  [Refe.  13,  14,  15,  &  16]  indicated  the  inadequacy  of  the  classical  laminate 
theory.  These  shear  deformation  theories  do  not  satisfy  the  conditions  of  zero  transverse 
shear  stresses  on  the  top  and  bottom  surfeces  of  the  plate,  and  they  require  a  shear 
correction  to  the  transverse  shear  stifGaess.  Three-dimensional  theories  of  laminates,  in 
which  each  layer  is  treated  as  a  homogeneous  anisotropic  medium  [Refe.  17,  18,  &  19] 
are  intractable  as  the  number  of  layers  becomes  moderately  large. 

Different  higher-order  laminated  plate  theories  have  been  proposed,  t\bich  account 
for  the  transverse  shear  strains.  Examples  of  such  theories  are,  Whitney  and  Pagano  [Ref 
20],  Whitney  et.  al  [Refe.  21,22],  Lo  et  al.  [Ref  23]  and  Nelson  et.  al.  [Ref  24]  In  these 
higher-order  theories,  an  additional  dependent  unknown  is  introduced  into  the  theory  with 
each  additional  power  of  the  thickness  coordinate. 

A  single  two-dimensional  theory  of  plates  that  accurately  describes  the  global 
behavior  of  laminated  plates  seems  to  be  a  compromise  between  accuracy  and  ease  of 
analysis.  A  single,  higher-order  theory  described  by  Reddy  [Ref  25],  is  such  a  theory, 
as  it  is  accounts  not  only  for  transverse  shear  strain  but  also  for  a  parabolic  variation  of 
the  transverse  shear  strains  through  the  thickness.  Consequently,  there  is  no  need  to  use 
shear  correction  coefficients  in  computing  the  shear  stresses.  This  theory  is  used  as  a 
prime  base  in  this  finite  element  analysis. 

The  finite  element  analysis  of  laminated  composite  plates  has  been  presented  by 
several  authors,  Reddy  et.  aL  [Refe.  26,36],  and  Noor  [Refe.37,38],  for  bending  and 
vibration  analyse.  Mawenya  [Ref  40],  developed  a  general  quadratic  isoparametric 
muhilayer  curved  plate  element,  and  Panda[Ref41],  presented  a  superparametric, 
quadratic  plate  element  for  the  plate  bending  analysis.  Fortier  and  Rossettos  [Ref  42], 
Sinha  and  Rath  [Ref  43],  analyzed  firee  vibration  and  buckling  of  thick  plates  of 
unsymmetric  cross-ply  construction.  Dong  [Ref  44]  has  given  the  solution  for  the 
dynamic  response  of  a  simply-supported  rectangular  plate. 

In  the  first  part  of  the  present  work,  a  finite  element  model  for  a  laminated  conposite 
plate  is  developed  using  a  simple,  higher-order,  shear  deformation  theory  with 
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Hamilton’s  principle  for  the  formulation  of  the  equations  of  motion  [Refe.  45,46].  A 
standard,  four  node,  rectangular  element  with  seven  degrees  of  freedom  at  each  node  is 
developed  for  the  analysis  of  a  flexible  laminated  plate  having  a  constant  thickness  for  any 
individual  layer.  The  displacement  model  is  so  chosen  because  it  can  represent  adequately 
the  parabolic  distribution  of  transverse  shear  stresses  and  the  non-linearity  of  in-plane 
displacements  across  the  thickness.  A  bending,  free  vibration  and  stress  analysis  problem 
is  examined  for  different  loadings,  boundary  conditions,  and  fiber  orientation  angles.  The 
results  are  coiiq>ared  with  existing  anafytical  and  numerical  solutions.  Hence  the  present 
element  formulation  denK)nstrates  its  applicability  for  a  wide  variety  of  laminate  conposite 
plates. 

Given  the  result  of  continuous  con:q)eting  requirements  for  in:q>roving  the  weight, 
interdisciplinary  performance,  tenperature  stability,  versatility,  and  reliability  of  propitlsion 
and  aeroqjace  conponents,  the  development  of  a  new  generation  of  composite  materials, 
called  ‘  intelligent/smart  conq)osite  materials  ‘  is  receiving  growing  attention,  which  is  the 
concern  of  this  dissertatioa 

The  reader  is  referred  to  books,  Cady  [Ref  47],  Dceda  [Ref  48],  and  Nye  [Ref  49], 
for  piezoelectricity  and  for  the  development  of  the  constitutive  relations  for  piezoelectric 
materials.  Tiersten  [Ref  50],  and  Rogacheva  [Ref  51],  contain  methods  for  solving  the 
differential  equations  of  the  theory  for  piezoelectrical  plates  and  shells,  respectively . 

Several  researchers  have  studied  the  interaction  between  the  mechanical  properties 
and  the  electric  field.  Crawley  et.  al.  [Ref.  52],  developed  piezoelectric  elements  for 
laminated  beams,  and  with  his  co-workers[Ref  53],  has  expanded  the  work  and 
considered  the  piezoelectric  actuators  to  be  plies  of  laminated  plate  and  used  the  Rayleigh 
Ritz  method  to  study  the  deformations  of  a  smart  plate.  They  also  modeled  the  induced 
strain  actuation  for  a  beam  [Ref  54],  and  for  a  truss  element[Ref  55].  Lee  [Refe.  56, 
57],  derived  a  theory  for  laminated  piezoelectric  plates  based  upon  classical  plate  theory. 
His  experimental  [Ref.  58],  results  showed  that  PVDF  or  PVDF2  actuators  can  generate 
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plate  bending  and  twisting  independently  or  simultaneously,  and  they  are  suitable  for 
active  d^^nping  control  of  a  flexible  structure. 

Tiersten  [Ref.  50]  modeled  single-layer  piezoelectric  plates,  including  the  charge 
equations,  but  did  not  study  laminates.  Wang  and  Rogers  [Ref  59],  appUed  the  classical 
laminated  plate  theory  to  model  laminated  plates  with  spatially  distributed  actuators. 
Tzou  et.  aL  [Ref  60],  derived  governing  equations  for  piezoelectric  shells  using  first  order 
theory.  Mitchell  and  Reddy  [Ref  61],  presented  a  hybrid  theory  for  enhancing  laminated 
plates  based  on  modeling  the  electric  potential  through  the  laminate  thickness  with  a  1-D 
finite  element.  Hagood  et.  al.  [Ref  62],  naodeled  the  effects  of  the  dynamic  coupling 
between  a  structure  and  an  electrical  network  through  the  piezoelectric  effect  for  a 
cantilever  beam.  State  space  models  were  developed  for  three  important  cases:  direct 
voltage  driven  electrodes,  direct  charge  driven  electrodes,  and  an  indirect  drive  case, 
where  the  piezoelectric  electrodes  were  connected  to  an  arbitrary  electrical  circuit  with 
voltage  and  current  sources.  Ray  et.  aL  [Ref.  63],  presented  the  exact  solutions 
for  a  two  dimensional  analysis  of  a  plate  composed  of  piezoelectric  material.  This  study 
led  to  exact  solutions  for  a  conqxisite  plate  with  piezoelectric  actuators  and  sensors  [Refe. 
64,  65].  Lee  [Ref  67],  developed  an  analytical  approach  via  state  space  equations. 
Heyliger  [Ref.  68],  developed  an  exact  solution  for  the  free  vibration  analysis. 

The  previous  solutions  do  not  provide  the  results  for  large  complicated  structures 
with  integrated  materials.  Thus,  the  necessity  for  approximate  techniques,  such  as  the 
finite  element  method  arises.  A  few  papers  have  been  developed  addressing  the  analysis  of 
intelligent  structures  by  the  finite  element  method.  Allik  and  Hughes  [Ref  69],  presented 
a  tetrahedral  finite  element  for  three-dimensional  electroelasticity.  Based  on  this  model, 
Tzou  [Ref.  70],  proposed  a  method  for  isotropic  plates  using  isoparametric  hexahedron 
solid  elements.  Kagawa  et.  al  [Ref  71],  developed  a  method  for  a  transversely  vibrating 
bar  with  electrostrictive  transducers.  McDearmon  [Ref  72]  and  Tzou  [Ref  73], 
develojied  a  method  for  single  cases  of  deformation.  Ha  et.  al.[Ref  74]  anal>^d  a 
composite  plate  by  using  a  three  dimensional  brick  element.  Both  elements  used  in  these 
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methods  made  the  problem  conqjlex,  costly  and  required  some  special  techniques  to 
overcome  inaccuracies  and  disadvantages  of  modeling  a  plate  with  3-D  elements.  These 
elements  display  excessive  shear  stifBiess  as  the  element  thickness  decreases.  The  two 
dimensional  quadrilateral  plate  element  developed  by  Hwang  et.  aL  [Ref.  75],  is  more 
efficient  than  solid  elements,  but  it  appears  to  have  restricted  modeling  capabilities. 
Hwang  et.  al.  [Ref.76],  proposed  a  model  based  on  classical  laminated  plate  theory,  which 
neglects  the  effects  of  transverse  shear  stresses,  and  hence  is  inadequate  for  the  analysis  of 
moderately  thick  composite  structures.  Chandrashekhara,  et.  aL  [Ref.  77],  developed  a 
model  based  on  the  first  order  shear  deformation  theory,  which  needed  a  shear  correction 
coefficient.  Static  analysis  or  an  intelligent  plate  was  presented  by  Ray  [Ref.  78],  using  a 
higher  order  shear  deformation  theory,  which  added  additional  dependent  unknowns  and 
made  the  problem  costly  to  solve. 

In  the  second  part  of  the  present  work,  a  finite  element  formulation  is  developed  for 
modeling  the  static  and  dynamic  response  of  laminated  conqiosite  plates  with  distributed 
piezoelectric  actuators  and  sensors  subjected  to  both  mechanical  and  electrical  loading.  A 
sinople  higher-order  shear  deformation  theory  with  Hamilton’s  principle  is  used  to 
formulate  the  equation  of  motion  of  the  system  in  which  the  piezoelectric  layer  is  treated 
as  a  normal  lamina.  The  displacement  model  represents  the  parabolic  distribution  of 
transverse  shear  stresses  and  the  non-linearity  of  in-plane  displacements  across  the 
thickness.  The  model  is  valid  for  both  segmented  and  continuous  piezoelectric  elements, 
wWch  can  be  either  surfece  bonded  or  embedded  in  the  laminated  plate.  The  piezoelectric 
layer  can  be  isotropic  or  orthotropic,  and  the  structure  core  can  be  anisotropic 
(graphite/epoxy,  etc.)  or  isotropic  (aluminum).  A  new,  four-noded,  bilinear,  isoparametric, 
rectangular  element  with  seven  mechanical  degrees  of  fireedom  and  one  electrical  degree 
of  fi-eedom  at  each  node  is  developed.  The  electric  potential  is  treated  as  a  generalized 
electric  coordinate  like  the  generalized  displacement  coordinates  at  the  mid-plane  of  the 
actuator  layers.  The  structural  deformations  due  to  electrical  and  mechanical  loads  are 
computed  and  compared  to  the  available  analytical  and  finite  element  resuhs. 
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The  third  part  of  this  dissertation  is  directed  toward  shape  control  and  optimization  of 
the  actuator  voltages.  Several  works  have  been  conpleted  in  vibration  control.  For 
exanple  Baz  et.  al.  [Ref  79],  Challopadhyay  [Ref  80],  Chandrashekhara  [Ref  77],  Park 
eL  at.  [Ref  76],  Tzou  [Ref  81,  82],  Anderson  and  Hagood  [Ref  83],  and  Hagood  et.  al. 
[Ref  84].  Wang  and  Fuller  [Ref  85,  86,  87,  &  88],  have  done  several  papers  for  active 
structural  acoustic  control  using  piezoelectric  actuators.  A  few  studies  have  been  done  in 
the  shape  control  and  optimizatioiL  Ghosh  [Ref  89]  showed  a  model  for  plate  shape 
control  by  using  PZT,  Agrawal  et.  aL  [Ref  90],  developed  a  mathematical  model  for 
deflection  using  a  finite  difference  technique  and  estimated  the  optimal  actuation  voltages. 
Kirby  [Ref  91],  and  Koconis  et.  al.  [Ref  92],  presented  an  anatytical  method  to  determine 
the  voltages  needed  to  achieve  a  specified  desired  shape  with  minimum  error  between  the 
actual  shape  and  the  desired  shape.  Based  on  Koconis  method,  Varadarajan  et.  al. 
[Ref  93],  showed  a  model  for  shape  control  of  a  laminated  plate. 

In  the  third  part  of  the  present  woric,  an  optimization  algorithm  based  on  the  finite 
element  technique  is  developed  to  determine  the  optimal  voltages  applied  to  the  actuator 
to  minimize  the  error  between  the  desired  shape  and  the  actual  shape.  The  error  flinction  is 
defined  as  the  mean  square  of  the  error  between  the  points  in  the  actual  sur&ce  and  the 
points  in  the  desired  sur&ce  for  each  element  of  the  finite  element  grid.  The  original  shape 
and  the  desired  shape  of  the  plate  are  specified.  Thus  the  objective  is  to  find  the  optimal 
actuator  voltage  to  get  minimum  error  between  the  desired  shape  and  the  actual  shape. 
The  analysis  is  based  on  small  deformation  theory,  therefore,  the  specified  desired  shape 
must  be  within  the  region  of  small  deformation  fi’om  the  original  shape. 

Two  Matlab  codes  were  developed.  The  first  code,  ‘COMPZ’  is  an  interactive  finite 
element  code.  It  is  able  to  solve  a  laminated  composite  plate,  with  or  without 
piezoelectric  layers,  subjected  to  mechanical  and  electric  loads  for  different  boundary 
conditions.  The  code  is  able  to  analyze  either  a  conq)lete  or  a  quarter  plate  analysis  to 
save  computational  time,  and  it  will  perform  a  bending,  fi'ee  vibration  analysis,  and 
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a  stress  analysis  for  the  individual  layers.  It  is  valid  for  different  types  of  material  for  both 
piezoelectric  or  laminated  layers  (  anisotropic,  isotropic,  etc.). 

The  second  Matlab  code,  ‘OPTSHP’,  is  able  to  analyze  the  structure  and  confute  the 
change  in  shape  due  to  mechanical  and  electrical  loads.  This  code  includes  a  loop 
containing  an  optimization  algorithm  coupled  with  a  Matlab  Optimization  Toolbox  to 
compute  the  voltage  applied  to  each  actuator  to  minimize  the  objective  function.  The 
numerical  results  are  presented  for  each  part  of  the  analysis. 
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n.  FINITE  ELEMENT  ANALYSIS  OF  A  COMPOSITE  PLATE 


A.  INTRODUCTION 


The  partial  differential  equation  governing  composite  laminates  of  arbitrary  geometry 
and  boundary  conditions  can  not  be  solved  in  closed  form  [Ref  95].  Analytical  solutions 
of  plate  theories  are  available  in  the  literature  (Le.,  Navier  solutions,  Levy  solutions).  The 
Rayleigh-Ritz  and  Galerkin  methods  can  also  be  used  to  determine  approximate  analytical 
solutions,  but  they  too  are  limited  to  sinple  geometries  because  of  the  difficulty  in 
constructing  approximation  fmctions  for  complicated  geometries.  The  use  of  numerical 
methods  &cilitates  the  solution  of  these  conplicated  problems.  Among  numerical 
methods,  the  finite  element  method  is  the  most  effective. 

There  are  several  types  of  finite  element  models  that  have  been  developed  for  plate 
theories.  These  can  be  grouped  into  three  major  categories;  displacement  models,  mixed 
and  hybrid  models,  and  equilibrium  models. 

The  displacement  finite  element  models  of  plate  theories  are  based  on  the  principles  of 
virtual  displacements,  where  all  governing  equations  are  ejpressed  in  terms  of  the 
displacements.  The  mixed  and  hybrid  finite  element  models  are  based  on  modified  or 
mixed  variations  of  the  plate  theories,  in  which  both  displacements  and  stresses  are 
independently  approximated.  The  equililsium  models  are  based  on  the  principle  of  virtual 
forces.  Among  the  three  types  of  models,  the  displacement  finite  element  models  are  the 
most  commonly  used  in  commercial  finite  element  programs. 

The  objective  of  this  chapter  is  to  develop  a  finite  element  model  for  a  laminated 
composite  plate,  using  a  single,  higher-order,  shear  deformation  theory.  The  model  is 
then  validated  for  different  boundary  conditions  and  loading  cases.  Bending,  fi'ee  vibration 
and  stress  analyses  are  performed  using  the  proposed  theory. 
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B.  GENERAL  THEROY 


Initial  configuration 
of  a  vertical  line 


Average  deformation 
configuration  of  the 
vertical  line 

Mid  surface 

Actual  deformation 
configuration  of  the 
vertical  line 


Figure  2.1 :  Deformation  of  the  normal  to  the  mid-plane  ‘Trom  Ref  [45].” 


A  typical  rectangular  laminated  plate  has  a  length  a,  width  b,  and  thickness  t.  The 
laminate  is  composed  of  a  number  of  perfectly  bonded  orthotropic  layers  (laminae),  which 
are  placed  sequentially  one  after  another.  A  coordinate  system  is  adopted  such  that  the  x-y 
plane  coincides  with  the  mid-plane  of  the  plate,  and  the  z  axis  is  perpendicular  to  the 
plane.  The  present  theory  is  formulated  based  on  the  following  displacement  fields  [Ref 
45]; 
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(2.1) 


U{x,y,z)  =  U{x,y,Q)  +  z^,{x,ySi)-^z^Cxix,yS>)  +  z^g,{x>yS>) 

V(x,y,z)  =  Vix,y,0)  +  z^/x,y,0)  +  zV,(x,  >^,0)  +  zV^(x,y,0) 
=  Fo+z^^  +  zV^  +  zV^o 


lF(x,:v.z)  =  FF(x,>’,0)  =  lFo 
where; 

C/q,  Fgand  IFo  are  the  displacements  of  a  point  on  the  reference  sur&ce  with  coordinates 
(xo,  yo,  Zo).  <f>^  and  are  the  average  rotations  about  the  Y  and  X  axes  respectively  of 

the  normal  to  the  mid-sur&ce  of  the  undeformed  plate,  and  z  is  the  distance  of  a  point 
from  the  mid-plane  along  the  Z  axis  as  shown  in  Figure  (2.1).  The  remaining  terms 
correspond  to  higher  order  rotations.  Applying  the  conditions  that  the  top  and  the 
bottom  surfrice  are  free  from  transverse  shear  stresses; 

r„(x,>;,±r/2)  =  0  and  r^(x,>^,±r/2)  =  0  (2.2) 

where  t  is  the  plate  thickness.  For  an  orthotropic  plate,  this  means  that  the  shear  strains 
are  zeros,  Le. 


€^{x,y;tt  12)^0  and  £^(x,>^,±r/2)  =  0;  (2.3) 

Differentiating  equation  (2.1)  and  substituting  in  equation  (2.3) 


13 


=  <^  /  ^  +  dfi^  /  dc 

+  3zV.(.  /  <*1.-1«!  =  0 

sj  =  *>,.  ±  <<;.  +  (3  /  4)(V^  +  ^  /  4t  =  0 

=  <^  /  <^ + ^  /  4^ 

^j®|z=±r/2  ~  ^  '*'  ^  ^lz=±//2  ~  ^ 

f>zlz=±r/2=  ±  +  (3  / 4)/V^  +  <^  /  <^  =  0 

By  solving  equations  (2.4)  and  (2.5),  we  get 
Co  =  0;  =  0; 

Sio  =  -^Uo  +  ^/^)  and  $-^0  =  -^(<^3-0  +  ^  /  4^) 
Substituting  equation  (2.6)  into  equation  (2.1),  we  get 

U (x,y,z)  =Uo+z  ^,0-  (^,o  +  ^  /  ^) 

V(x,y,z)  =  Fo  +  z  -  ^[7)  Uo  +  ^  /  4) 

W(x,y,z)  =  JV, 


(2.4) 


(2.5) 


(2.6) 


(2.7) 
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The  displacements  at  any  point  in  the  plate  are  given  in  terms  of  the  seven  unknown 
quantities  A,  and  dt/v! 

C.  STRAIN-DISPLACEMENT  RELATIONS 


The  linear  strain-displacement  relation  in  the  Cartesian  coordinate  system  is 


'dujac 

*  =  ^ 

{dUldi>^dVia:)  > 

By  substituting  equation  (2.7)  into  equation  (2.8),  we  get 


-  4(^wM^  + 

-  4(^w/4;^  +  /  3r^ 

*  zz  * 

-4(2<?^w/<3:4;  +  ^,o.j-  +  /3^^)  ' 

E„ 

-  4(7  /  tfid^ldK  + 

-  4(r  /  r)^(^/4^  +  (!i,o)  +  ^0/4; 

(2.9) 
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where; 


e>U,y, 

4  = 

^xy  ~  ^0  j'  ^  ^0^  > 

I 

4  =  ^xo+f^oy> 

^yz  ~  4o  '*’  ^0,y’ 

Kl^Ky 

f^y  ~  4o^> 

■^jy  ~  ^xO,y  ^y0,x’ 

Ki=^(^ldc  +  ^,,)lt^; 

Kl=-4i^w/^^+^^^,)l(3t^y, 

and 

4  =  -4(2^w/a^+4<,^  +  4,)/3/“);  (2.10) 

which  includes  linear  strains^  curvatures,  twists  and  other  higher  order  curvatures. 

D.  SRESS-STRAIN  RELATIONS 

For  anisotropic  lamina,  there  are  no  planes  of  symmetry  for  the  material  properties.  If 
there  is  one  plane  of  material  property  symmetry,  where  the  plane  of  symmetry  is  z=0  a 
material  is  termed  monoclinic,  and  has  13  independent  elastic  constants  .  For  a  monoclinic 
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The  corresponding  Zz  can  be 
This  gives  the  reduced  stress- 


^4  ' 

1 

II 

Q3Q3 . 
Qa  ’ 

^4,= 

=  Q.- 

Q3Q3 . 
c  ’ 

^33 

C^2 

=  Q2  — 

C43Q3 

c  ’ 

Ha 

c 

=  <^44- 

Q3Q3 . 

■  C33  ’ 

C55 

II 

C5. 

il 

Qs 

=  c  • 

'-'65’ 

and 

II 

(2.13) 


or  in  short  form; 

Q  ^  Q3 

for 

i,  j  =  1,2,4 

II 

for 

ij  =  5,6 

and  the  QcoeflScients  [Refe.  96, 97]; 


Qi  ~ 


■"13 


1 

E2EjA 


Cjj  - 


C33  - 


E2E2A 

£,£3/1  ’ 

t23,+U2,L>32 

^3  +  ^2l^23 

E2E2A 

E^E2A  ’ 

l-u,3t;3i 

£,£3^  ’ 

L>32  +  U,2Hi 

1223  +  t*2,Ui3 

E^E^A 

£,£2^  ’ 

1-0,2^21 

£,£221  ’ 

(2.14) 


18 


and; 


Qs  ~  ^13 » 

Qe  “  ^23  *» 

1  -  -  ^1^3  - 

£,£2^3 

wiiere  £,  ,£2  ,£3  are  the  Young’s  moduli  in  the  1 , 2,  and  3  directions,  respectively;  is 
the  Poisson’s  ratio  for  transverse  strain  in  the  j-direction  >^n  stressed  in  the  i-directioa 

For  orthotropic  lamina  there  »e  two  orthogonal  planes  of  material  property  symmetry 
for  the  material,  and  symmetry  will  exist  relative  to  the  third  mutually  orthogonal  plane. 
The  stress-strain  relations  equation  (2.1 1),  in  coordinates  aligned  with  principal  material 
directions  have  no  interaction  between  normal  stresses  and  shearing  strains;  Le., 
C,4  =  C24  =  C34  =  C4J  =  0,  Thus  the  stress-strain  relation  of  an  orthotropic  lamina  is; 

fc;,  c;2  o  o 

<^y  C^l  C22  0  0  0  Sy 

0  0  C44  0  0  (2.15) 

O’*  0  0  0  Cj5  0 

[0  0  0  0 

If  the  fibers  are  oriented  at  an  angle  0  with  the  x-axis  as  shown  in  Figure  2.2,  the 
transformed  stress-strain  relation  for  a  lamina  will  be 
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Figure  2.2:  An  off-axis  unidirectional  lamina. 


r  > 

O’, 

■Qh 

Qi2 

Q>4 

0 

0  ■ 

f  ^ 

Q2, 

Q22 

Q24 

0 

0 

►  — 

Q4. 

Q42 

Q44 

0 

0 

o-» 

0 

0 

0 

Q55 

Q56 

0 

0 

0 

Q65 

0.66. 

in  short  form; 

W]  =  \q\{£] 

where; 

Qii  =  Cji cos^  6  +  2(Cj2  +  2C44)cos^^sin^  6  +  Cjj  sin^  6 
Qi2  =  (Cii  +  ^22  -  AC^J^)coi6  sin^^+  Cj2(cos'‘^-i-  sin^  6) 


(2.16) 


(2.17) 


(2.18) 


Q,4  =  (C^i  -  2C44  -  Cj2)cos^0sin^+  (Cjj  -  +  2C44)cos^sin^^ 

Q22  =  C^iSin^^H-  2(0^2  +2C44)cos^ftm^^  +  C^cos^^ 

Q24  =  (^1  “  2C44  -  Ci2)cos^sin^^+  (0^2  -  Cjj  +  2C^^)oos^6&m6 
Q44  =  (Cii  +  “  2C^44)cos^  044(00$'’  +  sin^^) 

Qj5  =  Cjj  cos^  6  +  C^sin^^ 

Q54  =  (Cjj  -  C^)cosftin^ 

Qj4  =  Cj5sm^^+  Cjjcos^^ 


E.  STRESS  RESULTANT-STRAIN  RELATIONS 


Combining  equation  (2.9)  with  equation  (2.16)  and  integrating  layer-by-layer  over 
the  thickness,  the  following  relations  of  the  stress  resultant  are  obtained  [Ref  96]; 


=  £%,(l,z,2')dr;  (2.19) 

iNy,My,Py)=  [^^^tTy(l,Z,Z^)ct; 

(A^^,A4„P^)=  ^^cr^(l,z,z^)dz; 

iQ^,Ry,)^llay,i\,z^)dz. 


The  above  equations  can  be  set  in  the  matrix  form  as; 
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or  in  contact  form 


(w)  =  [d]{«") 


^ere ; 

[n]=Rigidity  matrix 

[a]  =Extensional  stiflSiess  matrix 

[5]  =  Coupling  stif&iess  matrix 

[D]  =  Bending  stif&iess  matrix 

[£l,[F],  and[H]  =  Higher  order  matrices 

which  are  given  for  n  layer  by: 


*=1  * 


(forj,i=l,...,6) 


(2.20) 


(2.21) 
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F.  VARIATIONAL  PRINCIPLE 


The  generalizd  fonn  of  Hamihon’s  principle,  vsiiich  states  that  the  variation  of  the 
Lagrangian  during  any  time  interval  must  be  zero,  has  the  form: 

J[^^<J(3  +  irVr  =  0,  (2.22) 

where  tj  and  t^  are  the  time  interval. 

The  Lagrangian  3  of  a  body  is  defined  by  the  summation  of  all  kinetic  energy  X  and 


potential  energy  ti . 

3=  [{%-K)dV 

(2.23) 

where; 

(2.24) 

*  = 

(2.25) 

Thus  the  Lagrangian  is 

3  =  iy  (2.26) 

\^ere; 

^  is  the  velocity  ( time  derivative  of  the  displacement  q); 

3  is  the  Lagrangian; 

V  is  the  volume. 

The  virtual  work  SfV  done  by  the  external  applied  force  is 

SfV  =  [[Sq]' {P,)ds  (2.27) 


where; 


s  is  the  surfece  areas  at  which  the  mechanical  load  is  applied 
P,  is  a  surfece  load  vector  (N/m^). 

Substituting  equations  (2.26)  and  (2.27)  into  the  Hamilton’s  eqviation  (2.22),  yields  the 
variational  equation  as: 


Since  all  the  variations  must  vanish  at  i  =  t,  and  /  =  /j,  by  substituting  equation  (2.17), 
and  taking  the  variation,  the  variational  equation  takes  the  form: 

lp{SgY{g}<iV*l{&Y[Qle)dV-l{s,Y{P,}ds  =  0  (2.29) 

G.  FINITE  ELEMENT  FORMULATION 

The  objective  is  to  define  the  degrees  of  fiwdom  Uo,Vo,<t>xo,<|>yo,Wo,Wx,  and  Wy  in  the 
plate  in  terms  of  nodal  displacements  and  rotations  by  using  a  bilinear,  isoparametric, 
rectangular  element  with  four  node.  Each  node  of  the  element  has  seven  degrees  of 
fi’eedom.  Similar  interpolation  fimctions  for  the  element  coordinates  x  and  y;  the  in-plane 
displacements  Uo  and  Vo  and  the  two  rotations  <|>xo  and  (((yo  were  used.  They  were  defined 

by; 

P  = 

i=l 

where; 

p  is  the  value  of  the  variable  at  any  point  in  the  element 
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(Xc-a/2,yc-b/2) 

(Xc-a/2,yo+b/2) 

(-1,*1) 

(-14) 

Actual  element 

Master  element 

^  =  2(x-x^)/fl,  Tj  =  2{y-y^)lb 
Figure  2.3:  Actual  and  master  element. 

Pi  is  its  value  at  node  point 

Ni  is  the  interpolation  function,  which  in  the  natural  coordinate  system  is; 

A^,  =  l/4(l  +  ^^,.Xl  +  7^)  (2.31) 

where; 

^  and  7  are  the  local  coordinates  of  the  point,  and  4  =  -1,14 -1  and 
7  =  “1,-144  for  /  =  1,..,4,  as  shown  in  Figure  (2.3) 

The  transverse  displacement  is  interpolated  using  a  non-conforming  shape  function  [Ref. 
98],  which  can  be  given  by; 
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w(x,:v) 


[/,  g,  K  fl  gl  K  h  g^  h  fA  gA 


(2.32) 


where  for  /  =  1  -►  4 

/  =  (1  / 8)  (1  +  ^4X1  +  +  m  (2-33) 

g,  =  (a/16)^  (l  +  - 1) 

h,  =  {b|\6)^  (1 + mf<y + - 1) 

where; 

^  =  2(x-xJ/a,  Tj-2(j-y^)/b,  (2.34) 

are  the  centroid  coordinates  of  the  plate ,  and  the  non-dimensional  coordinates  of 
the  nodes  are  (-1,-1),  (1,-1),  (1,1),  and  (-1,1),  respectively. 

The  nodal  displacement  vector  at  the  first  node  point  on  the  reference  surfece  is; 


“  [^01  IqI  ^x0\  ^01  ^.JTl  ^J'l] 

(2.35) 

and  the  element  displacement  vector  \q^  is; 

{9e}  =  [9i  Qi  Qi  ^aJ 

(2.36) 

The  generalized  displacement  vector  at  a  point  is 

{9)=[^0  K  </>x0  4>y0  %  ^.y] 

(2.37) 

Substituting  equations  (2.30)  and  (2.37)  into  equation  (2.10),  we  obtain  the  strain  vector 
at  the  mid-plane  in  a  matrix  form 
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The  first,  second,  and  higher  order  derivatives  of  the  interpolation  functions  \^ch  are 
expressed  in  equations  (2.38)  and  (2.39)  can  be  evaluated  with  respect  to  the  global 
coordinates.  An  additional  conqjutations  were  involved  ■winch  can  be  eiqiressed  as.  Let  ^ 
represent  one  of  the  interpolation  fiinctions  which  has  been  used  (Le.  ,  _/] ,  gj,  and/or 

h)- 

The  first  order  derivatives  with  respect  to  the  global  coordinates  are  related  to  those  with 
respect  to  the  focal  (or  element)  coordinates  according  to: 


-1 

di 

(2.40) 


where  the  Jacobian  matrix  [  J]  is  evaluated  using  the  ^proximation  of  the  geometry  wiiich 
may  take  the  form: 


or 


x=x,  +  ^a/2; 

y  =  yc  +  i7b/2 


x  = 

y  = 


X1+X4 

2 

yl•^y4 

2 


+  ^a/2; 
+  i7b/2 


(2.41) 


The  second  order  derivatives  of  Ti  with  respect  to  the  global  coordinates  (x,  y)  are 
given  by: 
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(2.42) 


f 

\^%] 

-N" 

i 

a^% 

af 

-[j=} 

.  ^ . 

A^. 

\ 

where; 


\dt])  \^) 

dj 

dqd^.  dr\d^ 


2— 
dr\  dri 

dj  df. 

drj  dt] 


w= 


ae 

^x 

£y 

arf^ 

^x 

dnd^ 

a^an 

(2.43) 


(2.44) 


Tte  elements  of  the  matrices  [j,]  and  [jj] 


are  computed  using  equation  (2.41). 


The  generalized  strain  at  a  point  is  related  to  the  reference  sur&ce  strain  as: 


where; 


(2.45) 
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(2.46) 


lOOOOzOO 
0  1  0  0  0  0  z  0 
0  0  1  0  0  0  0  z 
0  0  0  1  0  0  0  0 
0  0  0  0  1  0  0  0 


0 

0 

0 

z^ 

0 


0  z^  0 

0  0  z^ 

0  0  0 

0  0  0 

z^  0  0 


0 

0 

z^ 

0 

0 


The  reference  surfece  ( mid  plane )  strain  can  be  expressed  as; 

Thus; 

By  using  the  transformed  stress-strain  relation, 

io’]=[0]{4 

in  equation  (2.49) , 

W]-lQlBpU} 


(2.47) 


(2.48) 


(2.17) 


(2.49) 


H.  EQUATION  OF  MOTION 


By  substituting  equations  (2.48)  and  (2.49)  in  the  Hamliton’s  equation  (2.29)  gives 
the  ^^em  equations  of  motion: 


[{Sq.}’ ANY[NW{i.]*  l{Sq^\BY[B,"(\Q\B\B\dv{q,]- |{,%,} Vl'fn}*  =  0 

(2.50) 
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In  matrix  form,  the  equatioDS  of  motion  can  be  e^qrressed  as: 


+[<;]{,,)  =  {/v) 

>^i)ere  the  element  mass  matrix  is; 

K]= 


The  sh£q)e  function  matrix  [AT]  is  given  by  [Ref  47]; 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

o' 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

f> 

& 

h 

0 

0 

0 

K 

0 

0 

0 

I 


denotes  the  matrix  assemblage  for  the  element  four  nodes,  and 


(2.51) 


(2.52) 


(2.53) 


N  = 


/.  0 
0  /, 
0  0 
0  0 
0  0 
0  0 
0  0 


0 

0 

h 

0 

0 

0 

0 


0 

0 

0 

h 

0 

0 

0 


0  0  0 

0  0  0 

0  0  0 

0  0  0 

/,  0  0 

0  /j  0 

0  0  /j 


(2.54) 
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in  /j  and  I 2  flre  the  normal  and  rotational  inertia,  respectively. 


/.I  ' 


(2.55) 


\niiere  p‘  is  the  mass  density  of  the  i*  layer.  The  element  stif&iess  matrix  has  the  form; 


[«,]=  i[BY[B]'[QlB,\B]dV 
=  llBflDlB\dA 


(2.56) 


where  the  rigidity  matrix  [d]  for  n  number  oflayers  is  given  ly; 


[a]  =  Z  (for  j,  i  =1,. .  .6) 


(2.57) 


The  rigidity  matrix  D  is  given  in  equation  (2.20 ),  and  the  mechanical  excitation  force  is 
given  by; 


Thus  the  consistent  load  vector  at  node  i  is 

O' 


“  I  ^0 


0 

0 

0 

8i 

L*(, 


Ids 


(2.59) 


where;  is  the  intensity  of  the  load  per  unit  area. 
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Assembling  all  elemental  equations  gives  the  global  dynamic  system  equation: 


[w){«}+[a:]{?}  =  {/>„)  (2,60) 

I.  NUMERICAL  INTEGRATION 

The  elements  of  the  equation  of  motion  (2.S1),  which  are  giv»i  individually  in 
equations  (2.52),  (2.56)  and  (2.59)  can  be  integrated  numerically  using  Gauss 
Legendre  quadrature  (see  App.  A).  A  M  integration  technique  of  3x3  Gauss  points 
is  used  to  perform  the  integrations,  and  it  gives  satis&ctory  results  (comparing  to  the 
exact  and  available  finite  element  solutions)  for  all  plates. 


[A/,]=/^[ArnsrlAf]ii, 


i«i 


0 

0 

0 

0 

f> 


Si 


\J\d^Tj 


(2.61) 


(2.62) 


(2.63) 
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J.  FREE  VIBRATION  ANALYSIS 


A  free  vihratioii  analysis  can  be  performed  using  the  equation  of  motion  (2.60), 
wbere  the  right  hand  side  equals  zero  (the  work  done  by  external  forces  is  equal  to  zero 
for  the  free  vibration  analysis).  The  effects  of  the  degree  of  orthotropy,  of  different 
stacking  sequences,  aspect  ratio  (span-to-thickness  ratio)  of  the  structure,  numbers  of 
layers  and  orientation  angles  of  the  fibers,  on  the  natural  frequencies  and  mode  shapes 
are  demonstrated  in  the  validation  part. 

K.  STRESS  ANALYSIS 


Snoonwl  sK«»  ftfdisn 


Figure  2.4;  Local  discrete  least  squares  smoothing  “From  Ref.  [99].” 

By  solving  equation  (2.60)  for  a  static  deflection,  the  element  displacement  vector  is 
computed  ,  v^diich  can  be  substituted  into  the  generalized  strain  at  a  point  equation 
(2.48),  then  into  equation  (2.17)  to  get  the  stress  vector  at  any  point  in  the  structure. 
Stresses  are  obtained  at  2  x  2  Gauss  sampling  points,  as  shown  in  Figure  (2.4)  [Ref.  99]. 
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In  local  discrete  smoothing  [Ref.  99],  it  is  assumed  that  is  an  exact  least 

squares  fit  to  selected  values  of  at  the  Gaussian  integration  points  sho^vn  in 

Figure  (2.4).  (exact  least  squares  fit  is  one  ^\1uch  passes  through  all  the  sanpling  points 
and  is  therefore  interpolatory  in  nature).  Thus,  the  smoothed  comer  nodal  stresses 
CT,  ,CTj,  and  a^  may  be  obtained  fi-om  the  ejquession. 


1  ^ 

1  + 

1 

1 

2 

“2 

2 

2 

1 

1 

~2 

2 

"2 

2 

CTu 

1 

1 

2 

~2 

2 

2 

^iv 

2 

>4 

~2 

1 

T. 

(2.64) 


where  and  are  the  unsmoothed  stresses  at  the  Gauss  points 

(cTj  =  as  shown  in  Figure  (2.4).  The  smoothed  stress  values  are  then  modified 

by  finding  the  average  of  the  nodal  stresses  of  all  elements  meeting  at  a  common  node. 


L.  VALroATION: 


To  test  the  validity  of  the  finite  element  analysis  technique  and  to  establish  its  range  of 
£q>plicability,  numerical  examples  are  investigated. 

Exanq^le  1. 

A  three-  layer  (0°-90®-0°),  simply  supported,  square  plate  of  three  different  span  to 
thickness  ratios  (X.= 10,20, 100)  is  tested  .  The  total  thickness  of  the  0°  and  90®  layers  is 
the  same.  Layers  at  the  same  orientation  have  equal  thicknesses  (i.e.  ti=t3  and  ti+t3=t2). 
The  material  constants  are  as  follows; 
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Table  2.1:  Boundary  Conditions  for  various  different  conditions  (x  -  constant). 


Sin^fy  supported 

Clanped 

Free 

U^O; 

F  =  0; 

(7  =  0; 

F  =  0; 

C/9t0; 

II 

0 

II 

0 

j 

1^  =  0; 

*  0; 

w  =  0; 

35 

11 

p 

_ ! 

0; 

w^9t0; 

H'.y  =  0 

w^  =  0 

En  =  172.4  GPa  (25x10*  psi); 

E22  =  6.9GPa(10*psO; 

Gi2  =  Gi3=3.45  GPa  (0.5x10*  psi); 

G23  =  1.38  GPa  (0.2  X  10*  psi); 

V12  =  Vi3=0.25 

Aspect  Ratio(  length  to  thickness  ratio)=10, 20  and  100; 

The  plate  is  subjected  to  a  double  sinusoidal  load 
q=qosin(7tx/a)sin(jty/b) 

where  a*b=L;  the  plate  lengtL  Table  2.1  gives  the  boundary  conditions  on  the  edge 
X  =  constant.  Similar  statements  can  be  made  for  the  edge  ^ere  y  is  a  constant  by 
interchanging  the  subscripts  for  X  and  y.  A  4x4  mesh  quarter  plate  system  with  3x3 
sampling  points,  is  used  in  all  exaiEq)les.  The  results  are  given  in  Table  2.2  in  normalized 
quantities.  The  results  are  normalized  as  follows: 

=  (1  /  ; 

(aL,5>)  =  (1  / 

w  =  n^Qyi  /  12X*  tq^ ; 
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Table  2.2:  Three-Layers  cross  ply  (0°/90®/0®)  square  plate  [Ref.  45]  subjected  to  sinusoidal 
loading  (r,=r3  =  tIA,  1 12). 
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wiiere; 

S  is  the  reference  source  A  Present  finite  element. 

BFE2/PandaandNatarajan[Ref.  41].  C  FEl/Reddy  [Ref.  27] 

D  FE3/Mawenya  and  Davies  [Ref  40].  E  ISPQ/Moser  et.  aL  [Ref  104]. 

F  FE4/Phan  and  Red(fy  [Ref  33].  G  Exact/Pagano  and  Hatfield  [Ref  16]. 

CPT  Classical  plate  theory. 

and 

X  =  alt,  vsiiere  r  is  the  structure  thickness. 

Q  =  4Gj2  +  [£„  +  £22(1  +  2^12)] !  0  -  t^t>2,) ; 

w^=w(al2,a/2,0); 

CTj  =  /  2,a  1 2;tlf2) 

cTj  =  o>(a/2,fl/2^1/4); 

<^3  =  (0,04: 1/2); 

=  error  in  w  fi-om  exact  value 
=  error  in  tr,  fi’om  exact  value. 

E„  =  error  in  a,  firom  exact  value. 

E„  -  error  in  a,  fi’om  exact  value. 

Table  2.2  shows  that  the  numerical  values  of  all  quantities  are  converge  to  the  CPT 
results  by  increasing  the  span  to  thickness  ratio  X .  A  model  devebped  by  Reddy  [Ref 
27],  based  on  the  penahyAfNS,  (YNS  is  a  theory  of  Yang,  Norris  and  Stavsky)  with  eight 
node  quadrilateral  elements  and  five  degrees  of  fi«edom  per  node  is  perform  better  for 
small  span-to-thickness  ratios  X .  Ebments  developed  by  Panda,  Mawenya  and  Moser, 
either  have  too  many  structural  degrees  of  freedom,  or  have  unsatisfectory  performance 
conqjarison  to  the  present  model.  The  Phan  and  Reddy  element  gives  very  good  results 
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conqjared  to  the  exact  solution  for  the  thick  plate  (small  X ),  but  the  error  of  this  element 
increases  with  the  increase  the  span  to  thickness  ratio.  Thus  the  conclusion  is  that  for 
higher  span  to  thickness  ratios  one  can  use  the  present  element ,  and  for  lower  span  to 
thickness  ratios,  Phan’s  element  can  be  used  for  better  accuracy.  For  length  to  thickness 
ratios  equal  to  fifty.  Figures  (2,5),  (2.6),  and  (2.7)  show  the  different  stress  distributions 
for  a  composite  laminate  using  the  present  finite  element  method. 
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Normal  Stress  distribution  Of  Simply  Supported  Composite  Plate 


Figure  2.5:  Distribution  of  the  in-plane  normal  stress  Oxx  for  a  singly  supported 
laminated  square  plate  (0*’/90V0°)  subjected  to  double  sinusoidal  load. 
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Normal  Stress  distribution  Of  Simply  Supported  Composite  Plate 


Figure  2.6:  Distribution  of  the  in-plane  normal  stress  Cyy  for  a  singly  supported 
laminated  square  plate  (OVmVO®)  subjected  to  double  sinusoidal  load 
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Shear  Stress  distribution  Of  Simply  Supported  Composite  Plate 


Figure  2.7:  Distribution  of  the  transverse  ^ear  stress  Oxz  for  a  singly  supported 
laminated  square  plate  (0°/90®/0°)  subjected  to  double  sinusoidal  load. 
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Exanqjle  2. 


The  foUowing  material  constants  are  used  in  the  free  vibration  analysis  of  a  composite 
plate  [Ref.  46]; 


E,/E2  =  40;  Gi2=Gi3=0.6E2; 

G23  =0.5E2,  vi2  =  0.25 

The  dimensionless  fundamental  frequency  of  the  present  finite  element  is  con:^)ared 
with  different  methods  of  analysis  given  in  Table  2.3,  for  span  to  thickness  ratio  equal  to 
five.  The  effect  of  span  to  thickness  ratio  and  the  fiber  orientation  angles  of  the 
dimensionless  fundamental  frequency  are  given  in  Figure  2.8.  These  results  validated  the 
finite  element  code. 

Based  on  the  results,  the  fundamental  frequency  can  be  increased  by  increasing 
lamination  angle  up  to  45°  except  the  case  of  two  layers  in  which  it  decreases.  Increasing 
the  number  of  layers  with  fixed  thickness  increases  the  fundamental  fi^uency  due  to  the 
decreasing  of  flexural  extensional  coupling  effect.  The  dimensionless  frequency  increases 
as  long  as  the  span  to  thickness  ratio  decreases  and  reaches  a  maximum  value  at  a  fiber 
angle  equal  to  45°. 
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Table  2.3:  Dimensionless  fundamental  frequency  for  four  layers  (0°/90®/0°/90®)  with  span 
to  thickness  ration  equal  five 


Method 

Normalized  Natural  Frequencies’ 

Present  FE 

0.4500 

HSDT* 

0.44694 

HSDT" 

0.44686 

HSDT' 

0.44686 

FSDT" 

0.44083 

FSDT' 

0.45083 

Noor" 

0.42719 

__5 

0.66690 

CPT' 

0.66690 

Where; 

HSDT  is  the  hitter  order  shear  deformation  theory  [Ref  35]. 

FSDT  is  the  first  order  shear  deformation  theory  [Ref  35]. 

CPT  is  the  classical  laminate  plate  theory. 

*  Results  obtained  using  the  finite  element  solution  [Ref  35]. 

Results  obtained  using  the  Navier  Solution. 

‘  Results  obtained  with  the  exact  solmion  [Ref  35]. 

**  Results  obtained  by  applying  a  finite  difference  scheme  to  the  equations  of  the  three- 
dimensional  elasticity  theory. 

*  “From Ref  [35]” 
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Wn=w*sqrt(row*t''2/E2) 


Effects  Of  Fiber  Angles  &  Span  to  Thickness  Ratios  On  Dimensionless  Natural  Frequency 


Figure  2.8:  Effects  of  the  fiber  angles  &  span  to  thickness  ratio  on  the  dimensionless 
fundamental  natural  frequency  for  a  laminate  (0°/-6®/0°/-0°). 


m.  FINITE  ELEMENT  ANALYSIS  OF  A  SMART 
COMPOSITE  PLATE 


A.  INTRODUCTION 

The  rapid  development  of  high-speed  con5)uters  has  fecilitated  the  use  of 
computational  techniques  in  a  variety  of  engineering  applications.  The  finite  element 
method  is  one  of  the  most  popular  and  powerful  techniques  in  modem  engineering  design 
and  analysis  of  con^licated  structures  and  multifield  problems.  Most  of  the  research  on 
active  piezoelectric  stractures  has  focused  primarily  on  ejq)erimental  and  theoretical 
studies,  and  there  has  been  little  development  of  general  purpose  piezoelectric  finite 
element  codes  [Ref  99]. 

In  general,  experimental  models  are  limited  by  sizes,  cost,  and  other  laboratory 
unknowns.  Theoretical  models  can  be  more  general,  however,  analytical  solutions  are 
restricted  to  relatively  sin5)le  geometries  and  boundary  conditions.  In  the  case  of 
complicated  geometries  and/or  boxmdary  conditions,  both  theoretical  and  e)q)erimental 
techniques  encounter  technical  difficulties.  Thus,  finite  element  development  becomes 
very  inqjortant  in  the  modeling  and  analysis  of  a  elastic/piezoelectric  coupled  sy^em.  A 
typical  elastic/piezoelectric  stmcture  is  con^sed  of  a  main  elastic  continuum,  such  as 
aluminum  or  graphite  epojty,  with  coupled  (surfece  bounded),  or  embedded  piezoelectric 
sensors  and  actuators.  The  thickness  of  the  main  stmcture  can  be  about  two  or  three 
times  thicker  than  that  of  the  piezoelectric  layer.  Thus,  it  would  be  very  inefficient  and 
time-consuming  if  the  entire  stmcture  were  modeled  by  isoparametric  hexahedron  or 
tetrahedral  solid  elements.  Thus,  the  development  of  a  new  finite  element  model  applied 
to  piezoelectric  conqx)site  laminated  plates,  based  on  a  simple,  higher-order,  shear 
deformation  theory  is  essential  to  investigate  the  stmctural  response  due  to  the  applied 
field  and  ways  to  control  the  shape  of  the  stmcture. 
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In  this  cb^ter  a  finite  element  model  is  developed  and  a  single  higher-order  shear 
deformation  theory  derived  with  Hamilton’s  principle  used  to  formulate  the  equations  of 
motion  of  the  structure.  The  model  is  valid  for  a  piezoelectric  layer  either  surfece  bonded 
or  embedded  in  the  laminated,  plate.  The  piezo-lamina  could  be  patches  or  completely 
cover  the  sur&ce  as  shown  in  Figures  (3.1)  and  (3.2).  A  four  noded,  bilinear, 
isoparametric  rectangular  element  with  seven  mechanical  degrees  of  fireedom  and  one 
electrical  degree  of  fi'eedom  is  developed.  The  electric  potential  is  treated  as  a  generalized 
electric  coordinates,  similar  to  a  generalized  displacement  coordinates  at  the  mid-plane  of 
the  actuator  layer. 
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b.  Simply  supported  plate  with  piezoelectric  layer  cmnplete  cover  the  surfece 

Figure  3. 1  Composite  plate  with  piezoelectric  sensor  and 
actuator  conpletely  covering  the  surface 
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Piezoelectric  Actuators 


Piezoelectric  Sensors 


Figure  3.2:  Con^site  plate  with  distributed  segments  piezoelectric 
actuators  and  sensors 
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B.  VAWATIONAL  PRINCIPLE 


The  linear  piezoelectric  constitutive  equations  coupling  the  elastic  field  and  the  electric 
field  can  be  e3q)ressed  as  [Ref.  47]; 


_ 
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-^16 

^26 

^36 

^46 
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^66. 

.^16 

®26 

^36. 

(3.1) 


(3.2) 


Equation  (3.1)  describes  the  direct  piezoelectric  effect,  which  means  a  charge/vohage 
generated  by  an  imposed  force/pressure  to  a  piezoelectric  material.  The  converse 
piezoelectric  effect  is  described  by  equation  (3.2)  in  which  induced  stress/strain  are 
induced  due  to  an  externally  applied  vohage/charge.  The  equations  can  be  written  simply 
as: 


{D}=H^{4+MU}  (3.3) 

M=[c]{s}-[e]{E}  (3.4) 
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where; 

{ D)  =  Electric  displacement  vector  (Coulomb  /  meter  square) 

[e]  =  Dielectric  permittivity  matrix  (Coulomb  /  meter  square) 

(f)  =  Strain  vector 

[f"]  =  Dielectric  matrix  at  constant  mechanical  strain  (Farad  /  meter) 
(permittivity  conqwnent) 

[  E]  =  Electric  field  vector  (Volt  /  meter) 

{<t}  =  Stress  vector  (Newton/  meter  square ) 

{c}  =  Elasticity  matrix  for  a  constant  electric  field  (Newton  /  meter  square) 


It  is  assumed  that  the  principal  material  coordinates  coincide  with  the  coordinates  of 
the  problem  being  analyzed.  Thus  the  constitutive  relations  for  a  material  having 
orthohombic  mm2  symmetry,  including  piezoelectric  effects,  are  given  by: 
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In  this  analysis,  the  plane  stress  approximation  is  made  by  setting  <7^  =  0 ,  the  strain 
is  eliminated  from  the  constitutive  relations,  which  will  takes  the  following  form: 
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where; 

-  _  _  . 

®31  “  ®31  „  ®33> 

*'33 

—  ^23 

®32  ~  ®32  ”  ®33» 

*'33 

^24  ~  ^24’ 

^15  ~  ^15» 

^1  ~  ^1» 

^22  ~  ^2’ 
and 


—  %  + 

<^33 


(3.7) 


(3.8) 
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Table  3.1  Analogy  between  mechanical  and  electrical  quantities. 


Mechanical 

Electrical 

Force  density  Ps  (vector) 

Charge  density  p  (scalar) 

Displacement  q  (vector) 

Potential  O  (scalar) 

Stress  c  (second  -order  tensor) 

Flux  density  D  (vector) 

Strains  (second-order tensor) 

Electric  field  E  (vector) 

The  elastic  coefficients  Cjj  and  Cjj  can  be  obtained  in  the  same  way  as  Cjj  and  Cjj ,  which 

are  given  inequations  (2.13)  and  (2.14),  and  the  transformation  (if necessary )  can  be 
confuted  by  the  same  method  described  in  equation  (2.18).  The  analogy  between  the 
mechanical  and  electrical  quantities  [Ref.  69]  is  given  in  Table  3.1. 

The  Lagrangian,  3,  of  a  piezoelectric  body  is  defined  by  the  summation  of  all  kinetic 
energy  ,X ,  and  potential  energy  ,ti,  (including  strain  and  electrical  energies): 


:i=[{X-n)dV  (3.9) 

where; 

X  =  (310) 

and 

=  (3.11) 

Thus  the  Lagrangian  is 


(3.12) 


54 


w^ere; 


q  is  the  velocity  ( time  derivative  of  the  displacement  q); 

3  is  the  Lagrangian; 

V  is  the  piezoelectric  volume. 

The  virtual  work,  SfV ,  done  by  the  external  surfece  force  and  the  applied  surfece 
charge  density ,  // ,  ( CW)  applied  to  the  piezoelectric  body  is; 


(313) 

where; 

5,  and  $2  are  the  sur&ce  areas  at  which  the  mechanical  and  the  electrical  loads  are 
applied;  respective^ 

P,  is  a  sur&ce  load  vector  (N/m^) 

<P  is  the  electric  potential  (volt) 

The  minus  sign  in  eqmtion  (3.13)  occurs  because  in  the  variational  principle  for  the 
electromechanical  medium  it  turns  out  that  the  electric  enthalpy  H  =  U-  ;  takes  the 

place  of  the  internal  energy  function  i/  in  the  Lagrange  density;  Le.,  the  effective 
electrical  energy  content  of  /f  is  opposite  in  sign  to  that  of  U . 

By  using  Hamikon’s  principle 

=  (3.14) 

where  /,  to  tj  is  the  time  interval,  and  all  variations  must  vanish  at  t  =  r,  and  t  =  tj . 


55 


Substituting  equations  (3.12)  and  (3.13)  into  equation  (3.14)  yields  the  variational 
equation  as: 


By  substituting  equation  (3.3)  and  (3.4)  and  taking  the  variation  gives: 


(3.15) 


i: 


lU) + + {^]y]{E]w 


Wr  =  0 


(3.16) 


Since  all  the  variations  must  vanish  at  r  =  and  t  =  t2,  the  variational  equation  takes  the 
form: 


i\A^Y{i}  *  -  (&rw{£}  -  {«r^r{«}  -  {m'U']{E]yy 
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C.  FINITE  ELEMENT  FORMULATION 


To  generate  the  electro-elastic  matrix  relations  for  a  finite  element,  it  was  assumed 
that  the  surfiice  of  the  piezoelectric  layers  which  are  in  contact  with  the  laminated 
substructures  are  suitably  grounded.  Also,  since  the  thickness  of  the  piezoelectric  layers  is 
very  small,  it  is  reasonable  to  assume  that  the  electric  potential  fimctions,  \s1iich  yield  zero 
potential  at  the  interfece  between  the  actuator  and/or  the  laminated  substructure  and 
provides  a  linear  variation  across  the  thickness  of  the  sensor  or  actuator  layer  are  as 
follows: 

<t>\x,y,z)^{z-h^)<S>l{x,y)  (3.18) 

where  can  be  treated  as  the  generalized  electric  coordinate  similar  to  the  generalized 
displacement  coordinates  at  the  mid-plane  of  the  actuator  and  sensor  layers.  The 
generalized  electric  coordinates  at  any  point  within  the  element  can  then  be  e^ressed  in 
terms  of  i  nodal  variables  values  via  interpolation  fimction  N^\ 

♦.‘-[A'.K®;)  (3.19) 

where  is  the  nodal  generalized  electric  coordinate  vector  and  is  given  by 

{<«>.'}=[<!>.•  K  K  (3-20) 

with  0QI  ( i  =  1,..,4)  is  the  generalized  electric  coordinate  at  the  i*  node  of  the  element 

and  is  the  shape  fimction  matrix,  that  is,  the  same  linear  shape  fimction  used  for 
element  coordinates  x  and  y,  which  are  given  in  equation  (2.3 1). 

iV,=  1/4(1 +  #4)(1  +  7;x)  (3.21) 

where , 

=  N,  N,  W,]  (3.22) 
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Thus; 


(3.23) 


The  electric  field  vector  {£}  is  defined  by  the  electrical  potential  energy  by  using 
gradient  operator  V : 

{£}=-Vd>^  (3.24) 


Substituting  equation  (3.23)  in  equation  (3.24) ,  gives  the  electric  field  vector; 


where ; 


and 


[E]  = 


[s.]  =  v[w.] 


11 

_ 1 

^2 

^3 

^4 

0  0 

-iz-h^)  0 

0  -1 


(3.25) 

(3.26) 


(3.27) 


The  strain  vector  {f®}  at  the  reference  surfece  is  defined  by; 

(3-28) 

where ; 

is  the  element  displacement  vector  given  in  equation  (2.37). 
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[ B]  is  the  nodal  strain-displacement  matrix  given  in  equation  (2.38). 

The  generalized  strain  at  a  point  is  related  to  the  reference  surfece  strain  as; 


{^)=[5,F!  (3,29) 

Substitute  equation  (3.28)  into  equation  (3.29)  yields 

(3-30) 

where; 

lOOOOzOOO  Oz^O  o' 

OlOOOOzOO  0  Oz^O 
0  0  1  0  0  0  0  z  0  0  0  0  z^  (3.31) 

00010000  z^  0  0  0  0 

000010000  z^  00  0 

Substitution  of  equations  (3.25) ,  (3.30)  into  equation  (3.17)  yields; 


(  {ajlVn  />,}*,  +  ( {«>;)  V.fcr  -  =  0 


(3.32) 


D.  EQUATIONS  OF  MOTION 


1.  Plate  With  Actuators  Only 

The  variatiooal  equation  (3.32)  can  be  written  in  the  form; 

[AN\\N\dV%]{SiX  * 

(3.33) 

-  iKr[z,r!^‘KKK{®;K«'«}' 

- 1  [Af]' {?,}*,{<%,)'+  {«>'.)'  =  0 

Thus  the  equation  of  motion  in  matrix  form  is 

+[^w]ke}  =  {A/}  (3-34) 

[*,.,]{«,) +M{4-:)={«) 

After  adding  artificial,  linear,  viscous  damping  to  equation  (3.34),  the  equations  of  motion 
will  take  the  fom^ 

[^.]fe}  +  K,]fe}  +  [^,,]{9.}-[^^]{‘I>o}  =  M  (3.35) 
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where; 


\m.\  =  [AnXWv 

“  ljiNflm1iN]d(area) 


The  shape  function  matrix  [  A^]  is  given  by 


/=i 


and 


[m]  = 


N,  0  0  0 

0  iV;  0  0 

0  0  0 
0  0  0  JV; 

0  0  0  0 

0  0  0  0 

0  0  0  0 


/,  0  0  0  0 

0  /,  0  0  0 

0  0  0  0 
0  0  0  /j  0 

0  0  0  0  7, 

0  0  0  0  0 

0  0  0  0  0 
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1 - 
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0 

0 

0 

It 

J  i 

Si 

Si^ 

Si^ 

0  0 

0  0 

0  0 

0  0 

0  0 

h  0 

0  I2 


in  which  7,  and  I 2  are  the  normal  and  rotational  inertia,  respectively; 


(A,  A)  =  S 


where  p'  is  the  mass  density  of  the  /'*  piezoelectric  layer. 


(3.36) 


(3.37) 


(3.38) 


(3.39) 
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The  Hamping  matrix  |c^  J  is  defined  as  a  proportional  danping,  i.e. 


and 

(3.40) 

[*„]-iisnArnAj«]* 

- 

where  the  rigidity  matrix  [d]  for  n  number  of  layers  is  given  by 

(3.41) 

[^j]  =  S  f  '  Cy(l,zy  ,z^ ,z* ,z^)dz  (for  j,  i  =1,...6) 

and  its  elements  are  given  explicitly  in  equation  (2.20)  and 

(3.42) 

[*^]  = 

=  l\BY[Be2\B,)iA 

where  for  m  number  of  piezoelectric  layers  is  given  by; 

(3.43) 

*=i  ^  ■' 

and 

(3.44) 

[*-]  =  i[B^\[^,'{UY\B\BW  =  [*^]' 

(3.45) 
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and 
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To  specify  the  voltage,  the  distributed  applied  charge  density  required  to  create  a 
prescribed  electric  potential  distribution  on  the  surfece  of  the  actuator  layer  can  be 
determined  as; 

M(x,y)  =  -  (3.52) 

where  h^p and  hjpaie  the  distances  from  structure  middle  surfece  to  the  outer  and  the  inner 

surfrces  of  the  piezoelectic  actuator  layer,  respectively.  Assembling  all  the  equations 
gives  the  global  dynamic  system  equations 

[M]{4r)  +lcl{?)  +[a:„]|,)  -[v]f®)  = 

Note  that  the  mechanical  equation  is  coiq)led  with  the  electrical  equation;  in  which  {F}  is 
the  global  mechanical  excitation  and  {G}  is  the  electrical  excitation.  In  active  vibration 
control  application,  {G}  is  the  feedback  voltage  determined  by  the  control  algorithm. 
In  the  static  case,  equation  (3.53)  becomes; 

=  (3-54) 

K]W+[««]|®1  =  {G) 

By  performing  a  condensation  of  the  {d>}  degrees  of  freedom  from  equation  (3.54),  the 
equation  of  motion  will  be  as  follows: 

(3.55) 

where; 

(3.56) 
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and  the  electrical  potential  vector  is; 


(3.57) 


vrilich  can  be  used  to  calculate  the  voltage  distribution.  In  the  free  vibration  analysis,  { G} 
is  set  to  zero  so  that  the  voltage  distribution  associated  with  each  mode  can  be  estimated. 

Also  from  equation  (3.55),  the  feedback  force  jiy  |  can  be  expressed  as 

(3.58) 


2.  Plate  With  Actuators  And  Sensors 

For  actuator  and  sensor  the  variational  equation  (3.32)  will  takes  the  form: 

lAN'([N]dv{q,]{SqX 


65 


(3.59) 


Since  dq^  and  are  arbitrary,  equation  (3.61)  is  satisfied  only  if 


[t«,Ik}  +  [*Ml{®;L={«!  ‘Actuator’ 

After  aHHing  artificial  damping  and  assembling  all  the  equations,  the  global  dynamic 
system  equations  are; 

[M][q]  +[C]{^}  +[/:„]{?) =  {F}  (3.61) 

[A:c.,]j^}+[A^L{o}„  =  {G}  ‘Actuator’ 

+[^l{^}a  =  0  ‘Sensor’ 

where  {g} ,  {0}„,  and  {<!>},  are  the  global,  nodal,  generalized,  displacement  coordinates, 
and  the  global,  nodal,  generalized,  electric  coordinates  for  actuator  and  sensor, 
respectively. 
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The  global,  nodal,  generalized,  electric  coordinates  can  be  condensed  using  equation 
(3.61),  yielding  the  system  equations. 

(i/]{«j+[Cl{j)+[r]{,)  =  {f'}  (3.62) 

^\ilere; 

{r}  =  {F)  +  (3.63) 

Since  the  charge  applied  to  the  sensor  layer  is  zero  (  G=0),  the  voltage  from  the  sensor 
layer  can  be  written  as; 

(3-«5) 

The  values  of  the  matrices  in  equation  (3.60),  m^,  ,  k^^ ,  and  can  be 

conqjuted  from  equations  (3.36),  (3.41),  (3.43),  (3.45)  and  (3.46).  The  matrices  , 
and  can  be  computed  the  equations  (3.43),  (3.45),  and  (3.46)  for  the  sensor 
layers.  The  mechanical  and  electric  loads  are  confuted  using  equations  (3.49)  and  (3.51). 
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E.  NUMERICAL  INTEGRATION 


By  using  Gauss  Legendre  quadrature,  a  full  integration  technique  of  a  3x3  Gauss 
points  (see  Appendix  A)  are  used  to  con:q)ute  the  integration  of  each  element  of  the 
equations  of  motion  of  both  cases,  using  actuators  only  or  using  actuators  and  sensors 
(3.34)  and  (3.60).  Thus  equations  (3.36),  (3.41),  (3.43),  (3.45),  and  (3.46),  can  be 


nximerically  integrated  as  follows: 

l[NX[mlN\iA,  (3.66) 

(3.67) 

[t^]  =  I  [BY[BeZ\B^]u,  (3.68) 

[i«,]  =  f  [i!»]'[Z^'Z][s.]<i4,  (3.69) 
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0 
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0 
0 

W't  f,T>  ® 

f, 

g, 

.h. 

{g}  =  llN4u-h^)fdA,  (3.72) 

where  A*  is  the  master  element  area. 

F.  VALIDATION 

To  demonstrate  the  performance  of  the  finite  element  formulation  developed  in  the 
present  study,  the  numerical  results  were  compared  to  the  exact  solution  [Ref.  65],  and 
existing  finite  element  simulations  [Ref  78].  A  square  smart  plate,  consisting  of  a  three¬ 
layered  (0°/90“/0®)  cross-ply  laminated  plate  with  the  thickness  3  mm  was  used.  A  two 
piezoelectric  PVDF  layers  40  pm  each,  served  as  actuator  on  the  top  surfece,  and 
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a  second  as  sensor  on  the  bottom  surfece.  The  elastic  properties  that  simulate  a  high 
modulus  graphite/epoxy  composite  are  [Ref.  13]. 

Eii=  172.4  GPa  (25x10*  psi)  £22=6.9  GPa  (10*  psi); 

Gjj=Gj3=3.45  GPa  (0.5x10*  psi);  G23=1.38  GPa  (0.2x10*); 

Vi2^13=0.25 

The  piezoelectric  PVDF  layers  properties  are  [Ref  73]; 

Dielectric  permittivity 
e3i  =  0.0460 
e32  =  0.0460  Clm^ 
e33  =  0.0000  C/m'^^ 

Dielectricity 

el,  =  0.1062  X  lO'**  F/m 
4=0.1062x10*’  F/m 
4  =  0.1062  X  10*’  F/m 

Poisson  ratio  v  =  0.29 

Mass  density  p  =  0.1800  x  10^  kg/m  ^ 

Modulus  of  elasticity  E=2  x  10  ^  N/m  ^ 

The  mechanical  loading  and  the  electrical  potential  distribution  is  described  by: 

^  =  ^0  sin(mr  la)  sin(;^  /  b) 


and 


^ix,y,hop)  =  V  sin(m:  /  a)sin(;9;  /  b)\ 


where  is  the  intensity  of  load  per  unit  area  (N/m  \  and  F  is  the  anq)litude  of  O  in 
volts.  The  deflection  is  normalized  as  : 


w 


lOOEr 


w 


where; 

Ej-  is  the  transverse  Young's  modulxis  of  the  graphite/epoxy  layers 
A  is  the  span  to  the  thickness  ratio 
h  is  the  structure  thickness 
=10N/m'^^. 


Figures  (3.3),  and  (3.4  )  show  the  bending  deflections  verses  length  to  thickness  ratios 
for  a  smart  plate  subjected  to  a  double  sinusoidal  distribution  for  both  mechanical  and 
electrical  loads,  respectively.  Figure  (3.5)  shows  a  normalized  central  deflection  for  a 
sinqjly  supported  smart  plate  subjected  to  a  different  applied  voltage  values  A  three 
dimensional  plot  of  the  finite  element  grid  points  deflections  (81  nodes  )  under  the  same 
loads  is  shown  in  Figure  (3.6). 

From  the  results  we  can  conclude  that  a  new  finite  element  model  is  developed  using 
a  single  hjgher  order  shear  deformation  theory  to  analyze  a  smart  con^she  plate.  The 
numerical  results  are  conpared  with  the  exact  solution  [Ref  65],  and  existing  finite 
element  simulations  [Ref  78].  This  result  validates  the  present  finite  element  model.  The 
error  of  the  model  comes  in  to  the  picture  at  small  span  to  thickness  ratios  (3  and/or  4)  but 
it  decreases  dramatically  as  the  span  to  the  thickness  ratio  increases  and  it  attains  nearly 
the  exact  when  the  ratio  exceeds  one  hundred.  The  plate  deflections  can  be  increased  by 
increasing  the  applied  voltage  to  the  actuators. 
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Normalized  central  deflection 


Central  deflection  of  simply  supported  smart  composite  plate 


Length  to  thickness  ratio 

Exact/  [Ref.  65];  Reference  a,  FE/[Ref  78] 

Figure  3.3:  Bending  deflection  vs.  length  to  thickness  ratio  for  a  plate 
subjected  to  a  double  sinusoidal  mechanical  load. 
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Normalized  central  deflection 


Central  deflection  of  simply  supported  smart  composite  plate 


Length  to  thickness  ratio 


Exact/  [Ref.  65];  Reference  a,  FE/[Ref  78] 


Figure3.4:  Bending  deflection  vs.  length  to  thickness  ratio  for  a  plate 

subjected  to  double  sinusoidal  electrical  and  mechanical  loads. 
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Normalized  central  deflection 


Central  deflection 


Figure  3.5:  Normalized  central  deflection  vs.  applied  voltage  for  a  plate 
subjected  to  a  uniformly  distributed  electrical  load. 
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Deflection  W 


Grid  point  deflection  of  simply  supported  plate(PVDF, (0/90/0), PVDF) 


Figure  3.6:  Grid  point  deflection  for  simply  supported  plate  subjected  to 
a  double  sinusoidal  electrical  and  mechanical  loads. 
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VI.  SHAPE  CONTROL  AND  OPTIMIZATION 


A.  INTRODUCTION 


Future  technologies  ai  active  changes  in  the  structure  shape  such  as  airfoils 
spacecraft  antenna,  turbine  blades,  etc.  justify  and  necessitate  the  study  of  using 
piezoelectric  actuators  to  control  the  shape  of  structure  deformed  quasi-statically.  In  the 
previous  chapter,  a  finite  element  model  was  used  to  calculate  the  change  in  the  shape  of 
a  composite  plate  vdien  the  voltage  is  applied  to  the  piezoelectric  actuators. 

In  this  part  of  the  analysis,  a  new  model  is  developed  to  determine  voltage  applied  to 
the  piezoelectric  actuator  to  achieve  minimum  error  between  the  desired  shape  and  the 
actual  shape  with  prescribed  actuator(s)  placements.  The  mathematical  model  is  developed 
based  on  a  kinematical  relation  between  a  point  in  the  actual  sh^  and  a  corresponding 
point  in  the  desired  sh2q)e  for  each  element  of  the  finite  element  grid.  The  sj^em  of 
equations  for  the  error  function  is  formulated  based  on  a  finite  element  technique.  The 
optimization  algorithm  is  applied  to  the  error  (objective)  functions  through  a  finite  element 
analysis  program  called  OPTSHP.  The  code  performs  the  analysis  using  a  finite  element 
technique  and  determines  the  plate  deformation  due  to  the  ^plied  voltage  to  each 
actuator.  It  is  also  able  to  call  a  Matlab  optimization  function  to  confute  the  minimiun 
voltages  applied  to  each  actuator. 

B.  PROBLEM  STATMENT 

Consider  the  problem  of  a  plate  with  known  original  shape  and  actuator  configuration 
where  the  desired  shape  is  specified.  Thus  it  is  desired  to  find  actuator  voltages  to  achieve 
the  prescribed  shape,  which  minimizes  the  error  between  the  desired  shape  and  the  actual 
shape,  for  a  prescribed  actuator(s)  position,  to  achieve  a  minimum  error  (objective) 
function. 
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Piezoelectric  patch 


Figure  4,1  :  Zo  coordinate  shows  a  point  location  on  element  of  the  fimte 
element  grid 


The  analysis  is  based  on  small  deformation  theory.  Therefore  the  specified  desired 
shape  must  be  within  the  region  of  ‘  small  deformation’  firom  the  original  shape.  Also, 
the  analysis  uses  the  shape  change  of  a  reference  surfece.  Thus  within  the  approximations 
used  ,  the  shapes  of  either  the  ‘top’  or  the  ‘bottom’  surfece  of  the  element  can  be  taken 
as  the  reference  surfece. 

The  shape  of  each  element  of  the  finite  element  grid  is  described  via  a  suitably  chosen 
reference  surfiice.  The  shape  of  this  reference  surfece  is  defined  by  a  single  -coordinate 
of  every  point  element  sur&ce.  The  Zq -coordinate  is  perpendicular  to  the  x-y  plane  as 
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Actual  Sur&ce 


shown  in  Figure  (4.1).  The  shape  can  be  specified  a  polynomial  function  in  x  and  y 
such  as: 


^0=^0+  ^1^  +  +  ^3^  +  ^4^^  + 

(length  units) 

(4.1) 

Also  the  desired  shape  can  be  specified  by  the  same  way; 

( length  units) 

(4.2) 

^^ere  ai  and  bi  are  constants. 

C.  ERROR  FUNCTION  A 


Assume  that  a  given  set  of  voltages  is  applied  to  a  specified  number  of  actuators 
which  result  in  a  change  in  the  shape  of  the  plate  and  give  a  new  configuration  of  the 
structure,  (dashed  line  in  Figure  4.2.)  which  represent  the  actual  shape.  An  error  function 
will  be  introduced  which  includes  the  difference  between  the  actual  shape  and  the  desired 
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one  [Ref.  90]. 

The  error  function  will  take  the  form: 

A  =  '^dxdy  (length'’)  (4.3) 

where; 

A  is  proportional  to  the  distance  between  the  reference  surfece  of  the  desired  shape 
and  the  actual  reference  sur&ce 
Q,  is  the  element  domain. 

A  represents  the  difference  between  the  z-coordinate  on  the  reference  surfece  of  the 
desired  shape  and  the  z-coordinate  of  a  corre^nding  point  on  the  actual  reference 
surfece. 

The  error  function  can  also  be  the  sum  of  the  error  square  at  each  element 

A  =  (4.4) 

1=1 

where;  A  represents  the  difference  between  the  z-coordinate  of  a  point  on  the  desired 
shape  and  the  z-coordinate  of  a  corresponding  point  on  the  actual  surfece  for  each  element 
of  the  finite  element  grid. 

Thus  the  objective  is  to  make  the  error  function  A  as  small  as  possible  such  that  the 
desired  surfece  and  the  actual  surfece  are  nearly  identical  (Le.,yl  0)  for  a  prescribed 
actuator(s)  position  over  the  finite  element  grid,  with  optimal  applied  voltages  for  these 
prescribed  optimal  placements . 

From  Figure  4.3  the  distance  A  is  defined  as; 

A  =  (2o  +  &o)-(Zdes  +  ^des)  (length  units)  (4.5) 

To  evaluate  each  term  in  this  equation,  we  consider  a  nodal  point  A  on  the  original 
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Actual  Surfece 


Figure  4.3:  Difference  between  the  actual  surfece  and  the  desired  surfiice  in  the  x*z  plane 
‘TromRef.  [91].” 
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reference  surfece  at  (xo,yo).  Thus  Zo  is  the  z-coordinate  of  the  point  A  on  the  original 
reference  surfece,  is  the  difference  between  the  z-coordinate  of  point  A  on  the  original 
reference  surfece  and  the  z-coordinate  of  the  corresponding  point  A  in  the  actual  surfece  , 
which  can  be  ejqjressed  as : 


&o  = 


--^M  +  -^V  +  W 

^dx.  4^ 


(4.6) 


where ; 

«,v  and  w  are  the  displacements  of  the  node  A. 

w  is  the  displacement  tangential  to  the  reference  surfece  in  the  axis  x-z  plane; 

V  is  the  displacement  tangential  to  the  reference  surfece  in  the  off-axis  y-z  plane; 
w  is  the  deflection  normal  to  the  reference  surfece; 

z^  is  the  z-coordinate  of  point  B  on  the  desired  reference  surfece  at  (xo,yo). 


The  point  on  the  desired  reference  surfece  corresponding  to  point  A  on  the  actual 
reference  surfece  is  point  A  located  at  (Xq  +  ^o)-  ^  difference  between 

the  z-coordinates  of  points  A  and  B. 

It  be  approximated  by  the  first  two  terms  of  a  Taylor  series  expansion. 


= 


des 


V  3c 


des 


^  ) 


^0 


(4.7) 


The  distances  Sx^  and  4?^  th®  changes  in  the  x-  and  y-coordinates  between  points  A 
andB. 


(4.8) 
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<^0  = 


v-w- 


The  following  parameters  are  introduced  which  depend  on  the  geometry  of  the  original 
and  desired  sur&ce, 


7^=1  + 


( 


des 


4> 


(  ^o]  ( 


%= 


{  J 

By  using  this  parameter  the  objective  function  will  take  the  form: 

^  =  I,  [(^0  -  +  »£"  + 


(4.9) 


(4.10) 


D.  FINITE  ELEMENT  FORMULATION 

A  bilinear  isoparametric  rectangular  element  is  used  to  describe  the  shape  of  the 
element.  The  element  coordinates,  x  and  y  are  interpolated  using  a  linear  shape  function 
Nf ,  which  has  the  form: 

JV- 1/4(1  + ^4X1 +  »7^31)  (4.11) 


where; 

Ni  is  the  linear  interpolation  fimction. 
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^,7/  are  the  master  element  coordinates. 
x,y  are  the  actual  element  coordinates. 

The  actual  element  coordinates  x  and  y  can  be  transformed  to  the  master  element  by 
the  game  way  as  in  the  finite  element  analysis.  The  mathematical  relations  for  the  element 
transformation  will  have  the  form: 

x  =  Xc  +  4^l2;  (4.12) 

y  =  y,  +  Tjb/2; 
and 

^  =  2/o(x-xJ; 
rj  =  2lb(y-y,) 

By  substituting  equation  (4.12)  in  equations  (4.1)  and  (4.2),  the  actual  and  the  desired 
shape  equations  will  have  the  form; 

+  fl,(x,  +  ^a/2)  +  fl2(y^  + 176/2)  +  a^ix,  +  ^a/2)(y^  +  rjb/2)  (4.13) 
+  a^ix,  +  ^al2f  +  a^(y,  +  vb/2f 
and 

2des=bo+  ^(^c  +  ^«/2)  +  b2(yc  +  rjb/2)  +  b^ix^  +  ^a/2)(y,  +  776/2)  (4.14) 

+  64  (x,  +  ^al2f  +  b^iy,  +  776/2)^ 

(Zq  could  be  zero,  for  a  flat  plate,  or  a  first  order  polynomial ) 

By  the  chain  rule  of  partial  diflferentiation, 

3c  3^  &  3ri  3c 

&Q  _  ^ 

^  3rj  ^ 


84 


and 


S^des 

dc  dc  dt)  dc 

&des  _^des^  ^de,^ 

^  ^  dn  ^ 


(4.16) 


By  using  equations  (4.16)  and  (4.17),  the  three  geometric  parameters  qi,  ^2,  and  qs, 
which  are  e}q)ressed  in  equation  (4.9),  can  be  conq>uted.  The  objective  function  equation 
(4.10)  takes  the  form: 


-2<fe(^»7))  +  »ll(^.»7)w  +  »3!!(^.7)m  +  (417) 

The  electric  excitation  may  be  fectored  out  of  displacements  «,v,  and  w  in  eqiiation 
(4.17)  by  introducing  the  vectors  j7,vand  vv,  >^Wch  can  be  defined  as  follows.  The 
degrees  of  Ifreedom  of  the  system  can  be  determined  fi’om  the  equations  of  motions  (3.57) 
and  (3.64)  for  both  cases  of  actuator  and  actuator  and  sensor  that  are  given  in  chapter 
three  in  short  form  as: 

{?}=[r]''{f)  (4.18) 

where; 

{?}  isthearrayofthedegreesoffteedom 

[r]-'  is  the  inverse  of  the  structure  stiflfoess  matrix  (given  in  equations  (3.58)  and 
equation  (3.66)). 

{F*}  is  the  applied  electric  and  mechanical  loads  vector,  in  the  case  of  the  actuator 
only  which  is  described  in  chapter  three  as; 

{r)  =  {f)+[ji:^][/;;„,]'‘{G)  (4.i9) 
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in  case  of  applied  electric  load  only. 

{f"}  =[F^iF~.r(G) 

where  { G}  is  the  global  electric  actuating  force  in  which  for  one  element  g  is  defined  as; 

{«)  =  J[ 

and  the  distributed  applied  charge  density  is  described  as; 


(4.22) 


which  is  function  of  the  anqilitude  of  the  electric  potential  distribution  (F  in  volt). 

The  displacement  «,  v,and  w  are  defined  via  the  interpolation  shape  function  as: 
u  = 

/=1 

f=l 

4 

f=l 

«;,Vi,andWj  are  the  displacements  at  the  nodal  point  in  the  x,  y,  and  z  direction 
respectively,  which  can  be  picked  out  jfrom  the  deformation  array  {9}  expressed  in 
equation  (4.18)  such  as: 

«,=J^{F-}  (4.24) 

V,  =  ^{F') 

<*.  =  f,|f-} 
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where; 

is  the  flexibility  infliience  coefficient  rows  corresponding  to  the 

displacements  ^,v,,andH’^  in  the  global  matrix  inverse  [a^] 

By  substituting  u, ,  v,  and  w^  equation  (4.23)  takes  the  form: 

<4.25) 

<=1  /=1 

i=l 

/=i 

>^iere; 

fj  is  the  shape  function  used  for  the  deflection  w  at  node  j, 

Ni  is  the  shape  function  used  for  the  displacement  u ,  and  v  at  node  j,  and  w,v,and  w 

are  the  displacements  in  the  three  directions,  \^ch  are  confuted  fi'om  the  finite  element 
analysis  of  the  plate. 

These  displacements  are  a  function  of  the  voltage  ^plied  to  each  actuator.  To  fector 
out  the  voltage  in  the  objective  function,  u,v  and  w,  values  are  used  instead  of 
u,v,  and  w  in  equation  (4.17)  and  one  can  obtains; 

■^  =  |j  [(^o(^»'7) - ^do(^»7))  +  dxdy 

(4.26) 

which  is  a  function  of  the  an^litude  of  the  electric  potential  distribution  ( F  in  volt). 
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E.  NUMERICAL  INTEGRATION 

Equation  (4.26)  can  be  integrated  numerically  using  Gauss-Legendre  quadrature 
(see  Appendix  A).  A  full  integration  technique  of  a  3  x  3  Gauss  point  is  used  for  each 
element.  Thus,  the  error  function  for  the  one  element  can  be  written  as; 

(4.27) 

where  j  is  the  determinant  of  the  Jacobian  matrix. 

The  structure  objective  function  is  defined  as  the  summation  of  the  error  for  each 
element  such  as: 

i=l 

subjected  to 

Lower  limit  <  V  <  Upper  limit 

where  m  is  the  number  of  elements  of  the  finite  element  grid,  and  V  is  the  an^litude  of 
the  electric  potential  distribution  applied  at  each  actuator  (in  voltage). 

F.  OPTIMIZATION  ALGORITHM 

A  Matlab  function  /  min  and  /  mins  are  used  to  find  the  minimum  of  a  scalar 
function.  The  /  min  function  is  used  to  minimize  a  function  of  one  variable  on  a  fixed 
interval.  The  problem  is  mathematically  stated  as  : 
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(4.29) 


minimize  fix) 

X 

subjected  to: 

JC,  <  X  <  JCj 

where  fix)  and  x  are  scalars. 

The  /min  function  is  used  to  optimize  the  objective  function  in  the  a  case  of  constant 
{^lied  voltage  to  each  actuator  (F  is  scalar).  ' 

The  function  /  minj  is  used  to  find  the  mmimum  of  a  scalar  function  of  several 
variables  starting  firom  an  initial  estimate,  which  is  generally  referred  to  as  unconstrained 
nonlinear  optimization,  and  is  mathematically  stated  as: 

minimize  fix)  (4.30) 

X 

x^diere  x  is  a  matrix  or  vector  and  /  is  a  scalar  function. 

The  function  /  mins  is  used  when  dififerent  values  of  the  voltage  are  applied  to  the 
piezoelectric  actuators  (F  is  a  vector)  .  /  mins  uses  the  simplex  search  algorithm  of 
Nelder  and  Mead'°^,  for  minimization  of  a  function  of  n  variables,  which  depend  on  the 
comparison  of  function  values  at  the  (n+1)  vertices  of  a  general  simplex,  followed  by  the 
replacement  of  the  vertex  with  the  highest  value  by  another  point.  The  simplex  adapts 
itself  to  local  landscape  and  contracts  to  the  final  minimum.  The  method  is  shown  to  be 
effective  and  computationally  compact.  A  procedure  is  given  for  the  estimation  of  the 
Hessian  matrix  in  the  neighbourhood  of  the  minimtun  needed  in  statistical  estimation 
problems. 
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G.  SOLUTION  PROCEDURE 


The  flow  chart  of  the  solution  procedure  to  confute  the  optimal  voltages  applied  to 
the  actuators  and  to  minimize  the  objective  (error)  function  is  shown  in  Figure  4.4.  The 
procedure  is  as  follows: 

1.  Iiqjutthe  structural  data,  such  as  material  properties  for  both  graphite  epoxy 
and  piezoelectric  materials ,  plate  dimensions,  actuators  dimensions ,  boundary 
conditions  and  the  t^plied  voltages  to  the  actuator. 

2.  Assume  an  initial  guess  for  the  applied  voltages  to  each  actuator  at  prescribed 
selected  positions. 

3.  Call  the  optimization  tdgorithm  using  aMatlab  function  /minj  (in  case  of 
different  voltages  applied  to  different  actuators )  or  /  min  ( in  case  of  same  voltage 
applied  to  aU  actuators). 

4.  Compute  the  plate  deformation  due  to  the  current  applied  voltage,  using  the  finite 
element  analysis. 

5.  Evaluate  the  objective  fiinction  and  applied  voltages. 

6.  Check  the  convergence  of  the  objective  function  and  voltages.  If  they  converge, 
the  program  will  terminate. 

7.  Otherwise,  update  the  current  applied  voltages  to  the  actuator(s)  and  then  go  to 
step  three  and  repeat  the  procedure. 

H.  VALIDATION 

As  illustrated  in  the  flow  chart.  Figure  4.4,  the  optimal  location  of  the  actuator(s) 
position  to  get  a  minimum  error  function  at  a  specified  applied  voltage  is  determined. 
The  OPTSHP  code  for  the  optimization  algorithm  has  been  tested. 
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Figure  4.4:  Flow  chart  of  solution  procedure 
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Figure  4.5  Element  numbers  for  different  actuator  positions  for  a  sin:5)ly  supported 
square  con^site  plate  (a=  b). 


A  square  fiber-reinforced  con^sHe  plate  with  three  layers  (0**/90*’/0®)  and  one 
piezoelectric  layer  patch  on  the  top  surfiice  and  another  one  on  the  bottom  surfece  was 
used.  The  length  to  thickness  ratio  equal  fifty,  and  the  materials  properties  for 
graphite/epoxy  are  as  follows  [Ref  13]. 

Eii=  172.4  GPa  (25x10^  psi)  E22=6.9  GPa  (10^  psi); 

Gi2=Gi3=3.45  GPa  (0.5x10®  psi);  G23=1.38  GPa  (0.2x10®); 

Vi2=Vi3=0.25 


The  piezoelectric  PVDF  layers  properties  are  [Ref  73]: 
Dielectric  permittivity 

esi  =  0.0460  C/m''^ 
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632  =  0.0460  C/m^ 

633  =  0.0000  C/m^ 

Dfelectricity 

=0.1062  X  10-®  F/m 
4=0.1062x  10-’  F/m 
4=0.1062x10'’  F/m 

Poisson  ratio  v  =  0.29 

Mass  donsity  p  =  0.1800  x  10^  kg/m^^ 

Modulus  of  6lasticity  E=2  x  lO"^  N/m^^ 

Th6  mochanical  loading  is  takon  to  bo  zero,  and  tho  oloctrical  potontial  distribution  is 
assumed  to  be  uniform^  distributed  with  anqilitude  V .  The  deflection  is  normalized  as  : 

_  100£r 

H-  = - —W 

where; 

Ef  is  the  transverse  Young's  modulus  of  the  graphite/epoxy  layers 
X  is  the  span  to  the  thickness  ratio  (equal  fifty) 
h  is  the  structure  thickness 

A  nine  element  model  was  used  for  different  positions  of  the  actuators  as  shown  in 
Figure  4.5.  The  original  sh^  is  chosen  as  a  flat  plate  {z^  =  0)  and  the  desired  shape  is 

selected  as:  10“’ - 42  x  10“”j^  +  6 x  10"V^ • 
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Exan^le  1. 

In  this  case  the  error  between  the  actual  sh^  and  the  desired  shape  due  to  a  specified 
amount  of  the  applied  voltage  to  the  actuator  is  computed  by  placing  the  actuator(s)  at 
different  positions  on  the  finite  element  grid  (Figure  4.5  ).  The  error  function  values  are 
given  in  Table  4.1  for  an  ^plied  voltage  equal  to  100  volts  for  each  actuator  position. 

Exan^le  2. 

The  predicted  voltages  and  the  grid  points  transverse  deflections,  for  both  the  actual 
and  the  desired  transverse  deflections,  (in-plane  deformations  are  not  included)  are 
con?)uted  for  two  cases.  Figure  (4.6)  shows  the  first  case,  when  the  actuator  is  placed  at 
element  number  five.  The  second  case,  vriien  two  actuators  are  used  at  elements  two  and 
eight,  is  shown  in  Figure  4.7. 

Exan^le  3. 

The  optimal  value  of  voltage  applied  at  the  actuator(s)  to  obtain  minimum  error 
between  the  actual  shape  and  the  desired  shape  (minimum  error  fiinction)  can  be 
conq>uted  as  shown  in  the  fiow  chart.  Figure  4.4,  using  the  Matlab  f  min  function  with 

certain  upper  and  lower  limits  of  the  applied  voltage,  such  as  -100<  V  <  100.  The  optimal 
voltage  jqiplied  to  the  actuator  for  different  positions  on  the  grid  and  the  corresponding 
minimum  value  of  the  error  function  are  given  in  Table  4.2. 

The  reported  values  were  the  global  minimums.  This  was  acconqjlished  by  applying 
different  initial  values  of  voltages  (Le.  -500  to  500  Volt.)  to  the  actuator  and  determining 
the  numerical  value  of  the  objective  function  at  these  different  voltages.  The  minimum 
value  of  the  error  function  foimd  was  identical  to  the  value  computed  by  the  optimization 
algorithm. 
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Table  4. 1 :  Error  function  values  at  dififerent  actuator  positions  for  applied  voltage  equal 


100  vok. 


Actuator  Position 

Error  Function  Value 

Element  #  1 

f  =  5.87619804  e-18 

Element  #2 

f  =  5.05453907  e-18 

Element  #3 

f  =  4.44646994  e-18 

Element  #4 

f  =  IMOimS  e-18 

Element  #5 

f  =  9.30776749  e-18 

Element  #6 

f=  4.28947367  e-18 

Element  #7 

f  =  5.86952426  e-18 

Element  #8 

f  =  5.04220635  e-18 

Element  #9 

f  =  4.44063580  e-18 

•  The  applied  ventage  at  each  actuator  is  100  Volts. 


Example  4. 

In  the  case  of  four  actuators  used  at  a  time,  the  voltages  and  error  function  values  are 
given  in  Table  4.3.  Table  4.4  shows  the  results  for  two  actuators  used  at  a  time.  The 
results  in  both  cases  were  obtained  by  using  the  Matlab  function  /  min.v  with 
unconstrained  value  on  the  ^plied  voltage. 

An  optimization  algorithm  is  applied  to  get  minimum  applied  voltages  to  the  actuators  to 
minimize  the  error  function.  Tables  4.3  and  4.4  show  that  location  of  the  actuators  play  an 
in^rtant  role  in  minimizing  the  error. 
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Figure  4.7:  Actual  and  desired  transverse  deflection  at  the  grid  point 
for  actuators  at  elements  two  and  eight. 
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Table  4.2:  Optimal  iq)plied  voltage  and  error  function  at  different  actuator  positions. 


Actuator  Position 

Minimum  Applied 

Voltage  (Voh.) 

Minimum  Error 

Function 

Element  #  1 

V  =  78.226 

f=5.83495108e-18 

Element  #2 

V  =  70.944 

f=4.78995246e-18 

Element  #  3 

V  =  99.999 

f=4.44647064e-18 

Element  #4 

V  =  34.161 

f=5.88349351e-18 

Element  #  5 

V  =  37.346 

f=4.74690856e-18 

Element  #6 

V  =  75.062 

f=4.03165923e-18 

Element  #7 

V  =  78.611 

f=5.82972138e-18 

Element  #  8 

V  =  71.141 

f=  4.781 17803e- 18 

Element  #  9 

V  =  99.999 

f=4.44063651e-18 

Table  4.3:  Optimal  applied  voltages  and  error  fimctions  for  the  case  of  four  actuators  used 

at  a  time. 

Actuators  Positions 

Minimum  Applied 

Voltages  (Volt.) 

Minimum  Error 

Function 

Elements  #'  1, 3,  7  &  9 

Vi= -21.060 

Vs  =  98.104 

V7  = -16.473 

Vs  =  101.482 

f=  4.0531 138e- 18 

Elements  2, 4,  8  &  6 

V2  =  18.404 

V4  =  -29.429 

V6  =  65.262 

V8  =  22.349 

f=  4.0153 1283e-l  8 
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Table  4.4:  Optimal  applied  voltages  and  error  functions  for  the  case  of  pair  of  actuator 
used  at  a  time. 


Actuators  Positions 

Minimum  Applied 

Voltages  (Voh.) 

Minimum  Error 

Function 

Elements  #*  1  &  7 

Vi  =41.612 

V7  =  47.265 

f=  5.78007896e-18 

Elements  #’  2  &  8 

V2  =  35.943 

Vg  =  39.934 

f=4.74678213e-18 

Elements  #*  3  &  9 

V3  =  85.710 

Vs  =  90.547 

f=4.0583368e-18 

Elements  #*  1  &  3 

V,= -26.368 

Vs  =  175.055 

f=  4.101 14654e-18 

Elements  #*  4  &  6 

V4  = -8.231 

V6  =  79.660 

f=4.01544062e-18 

Elements  #*  7  &  9 

V,  =  26.0895 

Vs  =  175.193 

f=4.09288468e-18 

Elements  #*  1  &  9 

V,  =-17.007 

Vs  =170.053 

f=4.10970645e-18 

Elements  #*  3  &  7 

V3  =  169.712 

V7  = -16.138 

f=4.12087738e-18 
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V.  CONCLUSIONS 


A  finite  elen^nt  model  was  developed  to  analyze  the  behavior  of  a  smart  con^site 
plate.  This  intelligent  plate  was  con^sed  of  distribxited  sensor  and  actuator  layers  made 
of  PVDF  or  PZT  and  a  laminate  of  graphite/epoxy  layers.  In  the  analysis,  the  piezoelectric 
layer  was  treated  as  a  normal  lamina.  A  siipple  higher  order  shear  deformation  theory 
with  Hamilton's  princh)le  was  used  to  formulate  the  equations  of  motion.  A  four  node, 
bilinear,  isoparametric,  rectangular  element,  with  seven  degrees  of  freedom  at  each  node 
was  developed.  The  electric  potential  was  treated  as  a  generalized  electric  coordinate  like 
the  generalized  displacement  coordinates  at  the  mid-plane  of  the  actuator  and  sensor 
layers.  A  Matlab  code  ‘CMPZ’  was  developed  to  analyze  the  structure. 

Several  conclusions  can  be  made  from  the  results.  The  present  finite  element  method 
can  be  used  for  static  and  dynamic  analyses  for  both  thick  and  thin  con^she  structures 
with  distributed  piezoelectric  sensors  and  actuators.  The  numerical  results  generated  by 
the  developed  code  agree  very  well  with  the  exact  solution  and  other  published  finite 
element  solutions.  This  validates  the  finite  element  model  The  method  developed  is  much 
sinq)ler  to  formulate  and  more  efficient  than  models  based  on  hexahedral  or  tetrahedral 
solid  elements  and/or  three  dimensional  brick  elements.  The  error  of  the  method 
increased  for  small  span  to  thickness  ratios  and  decreased  dramatically  for  higher  span  to 
thickness  ratios.  A  Hermite  cubic  interpolation  function  was  used  to  approximate  the 
transverse  deflections,  however,  the  method  does  not  suffer  from  the  shear  correction, 
\^bich  is  problematic  in  the  first  order  shear  deformation  theory.  The  developed 
displacement  model  includes  the  parabolic  distribution  of  the  transverse  shear  stresses  and 
the  non-linearity  of  in-plane  displacements  across  the  thickness.  This  is  an  advantage  over 
the  classical  laminated  plate  theory,  which  neglects  the  effects  of  transverse  shear  stresses. 
The  number  of  degrees  of  freedom  of  the  element  used  in  the  present  method  is  one  third 
the  number  of  degrees  of  freedom  of  the  element  used  in  the  model  developed  by  the 
higher  order  shear  deformation  theory,  which  of  course  saves  con^utational  time. 
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For  shape  control  and  optimization,  an  optimization  algorithm  based  on  a  finite 
element  method  was  developed.  Kinematical  relations  were  used  to  formulate  the 
objective  (error)  function  as  the  summation  of  the  mean  square  error  between  a  point  in 
the  actual  sur&ce  and  a  corresponding  point  in  the  desired  surfece  integrated  for  each 
element  of  the  finite  element  gnd.  This  approach  is  an  in^rovement  over  the  method 
where  surfiice  error  is  determined  only  at  the  nodal  points. 

The  optimal  voltage  applied  at  each  actuator  to  achieve  a  specified  shape  with 
rriinimiim  eiTor  between  the  actual  shape  and  the  desired  shape  was  conqjuted  using  an 
optimization  code.  The  Matlab  code  ‘OPTSHP’  is  developed  in  conjunction  with  Matlab 
functions  which  gave  very  satisfectory  results.  The  code  is  able  to  determine  the  value  of 
the  error  function  for  different  actuator(s)  positions  for  a  specified  amount  of  the  applied 
voltage  to  the  actuator(s).  The  model  is  applied  for  both  cases,  either  constant  voltage 
applied  to  the  actuator(s)  or  different  voltages  applied  to  several  actuators.  The 
procedure  proposed  in  this  work  was  inqilemented  for  a  composite  plate  with  any 
demonstrated  number  of  piezoelectric  actuators.  It  was  demonstrated  firom  the  examples 
that  the  desired  shape  was  made  to  chosen  match  the  actual  sh^ie,  'wdiich  satisfies  small 
deformation  theory.  The  procedure  was  tested  for  more  general  layouts  of  actuators  and 
the  optimization  algorithm  was  applied  to  get  optimal  applied  voltages  to  the  actuators  to 
minimize  the  error  function. 

Future  work  is  proposed  in  the  following  areas; 

1.  Generalization  of  the  finite  element  model  to  non-rectangular  elements  and  shell 
elements. 

2.  Optimization  of  the  actuator  placements. 

3.  Closed  loop  shape  control  to  compensate  for  varying  external  mechanical  and  thermal 
loads. 
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APPENDEXA.  NUMERICAL  INTEGRATION 


For  the  isoparametric  element  considered  in  this  work,  an  exact  evaluation  of  the 
integration  for  mass,  stffiiess  and  the  equivalent  nodal  loads  was  not  always  possible 
because  of  the  algebraic  conqjlexity  of  the  differential  equations  representing  the  element 
used.  Since  the  interpolation  functions  are  easily  derivable  for  a  rectangular  element,  and 
it  is  easier  to  evaluate  integrals  over  rectangular  geometry,  we  transform  the  finite  element 
integral  statements,  defined  over  quadrilaterals,  to  a  rectangle.  In  such  a  case,  it  is 
necessary  to  use  some  numerical  integration  technique.  One  of  the  most  accurate  and 
convenient  methods  is  Gaussian  quadrature,  which  involves  approximation  of  the 
integrand  by  a  po^omial  of  sufficient  degree,  because  the  integral  of  a  polynomial  can  be 
evaluated  exactly  [Ref.  93].  If  an  integral; 

P=J['F(x)d&r  (Al) 

with  function  F(x)  is  approximated  by  a  polynomial; 

N 

(A.2) 

\\1iere  FJ  denotes  the  value  of  F(x)  at  the  I*  point  of  the  interval  (xi,X2)  and  are 
polynomials  of  degreeA-1.  The  representation  can  be  viewed  as  the  finite  element 
interpolation  of  F(x) ,  where  F]  is  the  value  of  the  function  at  the  I*  node.  For  our 
rectangular  element  the  integral  will  take  the  form; 

M  N 

[  F(x,y)dx<iy=l  ^ 

"  '  /=1  y=l 

(A3) 
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where; 

is  the  actual  element  domain  with  its  x,  and  y  coordinates  (  Fig.  2.3) 
is  a  linear  master  element  with  its  r  and  s  coordinates  (Fig.2.3) 

(^,  rf)  denote  the  Gauss  points  in  the  ^  and  17  directions.  ( Fig  2.3) 

R,  and  Rj  denote  the  corresponding  Gauss  weights  ( see  Table.  A.1 ) 

A/ and  N  denote  the  number  of  quadrature  points.  ( in  most,  cases,  the  interpolation 
functions  are  of  the  same  degree  in  both  %  and  ti  direction,  Le.  M  —  N . 

I  J1  is  the  determinant  of  the  Jacobian  matrix  J,  in  ^^ch  ,|  >  0  every\^ere  in  the 

element  ^2^. 

Geometrically,  the  Jacobian  represents  the  ratio  of  an  area  element  in  the  actual 
element  to  the  corresponding  area  element  in  the  master  element: 

dA  =  dxcfy  =  \j\d^7j  (A.4) 

The  number  of  Gauss  points  is  based  on  the  formula; 

N  =  mteger(lf2(p  + 1))  (A.5) 

which  means  that  the  smallest  integer  greater  than  1/2  (p  + 1) ;  >^iiere  p  is  the  degree  of  the 
polynomial  which  may  be  integrated  exactly  enqjloying  equation  (A.5). 

Table  A.1  gives  the  integration  order,  weighting  fector  and  the  location  of  the  Gauss 
points  for  linear,  quadratic,  and  cubic  elements  [Ref  94].  The  maximum  degree  of  the 
polynomial  refers  to  the  degree  of  the  highest  polynomial  in  4  and  t)  that  is  present  in  the 
integral  F{^,Tf)  of  equation  (A.3). 
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Table  A.1: 

Weights  and  Gauss  points  for  the  Gauss-  Legendre  Quadrature. 


N 

1 

1 

0.0 

2.0 

2 

0.5773502692 

1.0 

3 

0.7745966692 

0.5555555556 

0.0 

0.8888888889 

4 

0.8611363116 

0.3478548451 

0.3399810436 

0.6521451549 

5 

0.9061798459 

0.2369268851 

0.5384693101 

0.4786286705 

0.0 

0.5688888889 

6 

0.9324695142 

0.1713244924 

0.6612093865 

0.3607615730 

0.2386191861 

0.4679139346 

7 

0.9491079123 

0.1294849662 

0.7415311856 

0.2797053915 

0.4058451514 

0.3818300505 

0.0 

0.4179591837 

8 

0.9602898565 

0.1012285363 

0.7966664774 

0.2223810345 

0.5255324099 

0.3137066459 

0.1834346425 

0.3626837834 

N=2,  Linear,  N=3,  Quadratic;  N»4  Cubic,  etc. 

N  X  N  =the  order  of  integration 
%  -  Location  of  mtegraticm  points  in  master  element. 
R*  Weighting  fector. 
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APPENDEXB.  COMPUTER  PROGRAMS 


Two  Matlab  codes  were  developed,  'COMPZ'  and  'OPTSHP'.  The  'COMPZ'  code 
is  a  finite  element  analysis  program  for  a  smart  conqtosite  plate.  The  program  is  able  to 
compute  the  Static  and  the  dynamic  response  of  the  structure,  bending  deformation,  free 
vibration,  and  stress  analysis  subjected  to  both  mechanical  and  electric  load  for  various 
boimdary  conditions.  The  program  is  able  to  solve  a  composite  plate  with  and  without 
piezoelectric  layers.  The  composite  laminates  could  be  graphite  epoxy  with  reinforced 
fiber  or  ahimimim  layers.  The  piezoelectric  elements  can  be  segmented  or  continuous 
elements  and  can  be  either  surfoce  bonded  or  embedded  in  the  laminated  plate.  The 
procedure  to  run  the  program  is  as  follows: 

1.  Mesh  generation : 

Input  the  number  of  point  in  x  direction  Nx 

Input  the  number  of  point  in  y  direction  Ny 

Input  the  plate  length  in  X  direction  Lx 

Ii^ut  the  plate  length  in  y  direction  Ly 

Input  the  plate  thickness  in  z  direction  Lthk 

Input  the  number  of  patches  their  lengths  in  x  direction  (if  any) 

Input  the  number  of  patches  their  lengths  in  y  direction  (if  any) 

The  program  will  generate  the  mesh  automatically  and  compute  the  elements  numbers, 
nodes  coordinates,  and  the  elements  conductivity. 

2.  Structures  Data: 

Input  type  of  the  structures,  singly  supported  or  cantilever  plate. 

Input  materials  constants  for  the  structures  core  (graphite  epoxy  or  aluminum )  and  for 
the  piezoelectric  material.  E),  E:,  etc. 
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3.  Applied  Loads: 

Input  the  intensity  of  the  load  per  unite  area  q©. 

Ii^)ut  the  amplitude  of  the  {q)plied  voltage  V. 

4.  Method  of  analysis: 

Input  quarter  plate  analysis. 

Input  full  plate  analysis. 

Con:q>osite  plate  analysis. 

Smart  con^site  plate  anafysis,  with  piezoelectric  materials. 
Complete  cover  the  sur&ce,  or 
Uniform  patches  cover  the  sur&ce 

5.  Resuh 

Plate  deflection 

Natural  frequency  and  mode  shape 
Stress  value  for  each  layer 


The  "OPTSHP"  code  is  an  optimization  program  using  the  finite  element  technique  in 
which  the  optimization  algorithm  is  applied  to  the  objective  function  by  using  a  Matlab 
function  /  mins  through  finite  element  analysis  subroutines.  The  code  is  able  to  compute 
the  plate  deformation  due  to  the  applied  voltage  and/or  mechanical  load  and  determine  the 
optimal  placement  of  actuator.  Also  the  code  is  able  to  compute  the  minimum  amoimt  of 
the  applied  voltages  to  each  actuator  by  using  a  Matlab  function  /mins  and  /min. 

The  steps  to  run  the  code  are: 

1 .  Mesh  generation : 
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Input  the  number  of  point  in  x  direction  Nx 
Ii^)ut  the  number  of  point  in  y  direction  Ny 
li^ut  the  plate  length  in  X  direction  Lx 
Input  the  plate  length  in  y  direction  Ly 
li^ut  the  plate  thickness  in  z  direction  Lthk 
Input  the  number  of  patches 
Input  the  patches  element  numbers 

2.  Structures  Data: 

It5)ut  type  of  the  structmes,  singly  supported  or  cantilever  plate. 

Input  materials  constants  for  the  structures  core  (graphite  epoxy  or  aluminum )  and  for 
the  piezoelectric  material.  Ei,  E2,  etc. 

3.  Applied  Loads: 

Input  the  intensity  of  the  load  per  unite  area  qo  (if  any) 

Ii^ut  the  initial  value  of  the  applied  voltage  V . 

4.  Method  of  analysis: 

Input  full  plate  analysis. 

Smart  composite  plate  analysis,  with  piezoelectric  materials, 

Uniform  patches  cover  the  surfece 

5.  Determine  the  patch  location: 

Input  initial  guess  for  patch  position  and  con:q)ute  the  objective  function 

The  program  will  do  finite  element  analysis  and  compute  the  plate  deflection  due 

to  the  applied  voltage  at  this  place. 

Input  the  second  guess  at  other  location  of  the  finite  element  grid. 
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Repeat  this  dieck  and  corapaie  the  values  of  the  objective  function  for  each  case  and 
pick  the  minimum  value  which  represent  the  optimal  place  to  get  minimum  error 
between  the  actual  and  the  desired  shape. 

6.  Determine  the  optimal  amount  voltages  {q>plied  to  each  actuator: 

Input  initial  value  of  the  applied  voltages  for  the  each  actuator  at  the  selected  position. 
Call  the  optimization  algorithm  using  the  Matlab  function  f  mins  and 
/  min ,  they  will  perform  a  finite  element  analysis  for  this  value  of  the  voltages. 

The  code  will  update  the  voltage  values  and  repeat  the  analysis  again. 

The  code  will  check  the  global  minimum  of  the  objective  function  and  will  stopped  at 
the  global  minimum  with  minimum  amount  of  voltages  s^plied  to  each  actuators. 
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